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GODUNOV METHOD AND COMPUTER PROGRAM TO DETERMINE THE 
PRESSURE AND FLOW FIELD ASSOCIATED WITH A SONIC BOOM FOCUS 

By Lee W. Parker and Robert G. Zalosh 

Mt. Auburn Research Associates, Inc. 

Newton, Massachusetts 

SUMMARY 

A numerical method has been developed to calculate the 

flow field associated with sonic boom focusing. The computer code 

used in the calculations is capable of following weak N-waves (rela- 

-3 

tive overpressures of the order of 10 before focusing) for large 
distances as they approach the focus, in addition to providing the 
flow field at the focus itself. 

Results are presented for two types of problems. In one 
type of problem, the refraction and subsequent focusing of the N-wave 
is caused by a localized cold-spot in the atmosphere. In the second 
type of problem, the N-wave is assumed to be initially refracted (for 
example, by atmospheric inhomogeneities or aircraft maneuvers) into a 
prescribed concave shape. 

Several sample problems in each category were run. Result- 
ing overpressures at the foci range from 4.4 to 20 times the nominal 
overpressure. The typical length scale of the high-pressure focal 
region is of the order of one wavelength. These results are for 
hypothetical situations, not necessarily typical of supersonic air- 
craft booms. However, the computer code is now available for use 
with data taken from specific maneuvers and/or atmospheric disturbances. 

An interesting result of this investigation is the resolu- 
tion of the controversy concerning wave folding at a focus. The 
theory of geometric acoustics predicts that a concave shock front 



will fold over upon itself as it propagates through a focus (ref. 2). 
As opposed to this, Whitham (ref. 9) has claimed that a concave shock 
will straighten out without folding over. It appears from the calcu- 
lations reported here that both phenomena occur, but under different 
conditions. A weak shock wave with a relative overpressure much less 
than unity folds over, whereas a strong shock with a relative over- 
pressure of the order of unity (or higher) tends to straighten out. 



1. INTRODUCTION 


The overpressures in a sonic boom N-wave can be intensi- 
fied through the focusing phenomenon associated with concave shock 
fronts. Such concave shock fronts may be produced by atmospheric 
inhomogeneities and by aircraft maneuvers. This report is addressed 
to the computation of the flow fields that occur during focusing. Of 
particular interest is the calculation of the pressure at the focus 
as well as the extent of the high-pressure region surrounding the 
focus . 

The far-field disturbance from a supersonic aircraft is 

in the form of an N-wave with two shocks, "fore" and "aft," separated 

by a linear rarefaction wave. Since the overpressures in a sonic boom 
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N-wave are typically of the order of 10 times atmospheric pressure 
(i.e. extremely weak) before focusing, the theory of geometric acoustics 
is often employed to describe the approach to the focus. However, as 
explained later, geometric acoustics must be abandoned at the focus 
itself. On the other hand, standard numerical finite-difference 
techniques cannot be used to follow weak shocks for large distances 
because their artificial viscosity tends to smear out the shocks. 

The shock-following method developed during this study combines the 
advantages of both geometric acoustics and hydrodynamics. It preserves 
the weak shocks of the N-waves throughout the entire flow field, and 
accurately computes the pressure in. the focus. 

Throughout this report, both line foci and point foci will 
be discussed. Axisymmetric (r,z) geometry leads to a point focus, 
whereas two-dimensional (x,y) geometry leads to a line focus. For 
simplicity, the discussion to follow treats single shocks. It is 
understood that similar reasoning may be applied to both shocks of 
the N-wave. 

The first mechanism for generating concave shock fronts 
that will be discussed here is the refraction of a shock wave through 
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a localized cold region in the atmosphere. Consider the situation 
illustrated in Figure 1-1, where the bow shock from a supersonic air- 
craft passes through a spherical atmospheric "cold spot" (leading to 
a point focus as described below). If the sound speed variation with 
altitude is neglected, the shock wave propagates at a uniform speed 
before it penetrates the cold spot. Since sound speed in the cold 
spot is lower than the ambient sound speed, the portion of the shock 
passing through the cold spot is slowed relative to the rest of the 
shock. The resulting concave shock front focuses in such a manner 
that a shock cusp is produced at the focus. Thus, the action of the 
cold spot in producing a focus is similar to that of an optical lens. 

As Figure 1-1 indicates, beyond the focus, the depth of the cusped 
region decays asymptotically with time. 

The theory ordinarily used to predict sonic boom propaga- 
tion (ref. 1) is based upon a modified form of geometric acoustics. 
According to geometric acoustics, the wave front propagates along 
rays which are everywhere perpendicular to the wave front. A geo- 
metric acoustics description of the propagation of a concave shock 
is illustrated in Figure 1-2, where the rays are drawn, and Figure 1-3, 
where successive shock fronts are drawn. The caustic sheets shown in 
Figures 1-2 and 1-3 are defined as imaginary surfaces along which 
adjacent rays cross. The two caustic sheets shown in Figure 1-3 
correspond to a two-dimensional (x,y) geometry in which the focus is 
a line focus. In axi symmetri c (r,z) geometry, the caustic sheets 
would become a single surface of revolution, and the focus would 
become a point focus. 

One of the premises of geometric acoustics is that the 
local amplitude of the wave front is inversely proportional to the 
square root of the area of the ray tube formed by adjacent rays. Con- 
sequently, geometric acoustics predicts infinite overpressures at a 
caustic, where the ray tube areas vanish. Of course, non-linear dif- 
fraction effects, which are not accounted for in geometric acoustics. 
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limit the overpressure amplification to some finite value. This value 
may be computed by the method of the present report. 

Although geometric acoustics fails to predict shock over- 
pressures in the immediate vicinity of a caustic sheet, there is no 
reason to discount the qualitative geometric acoustics description in 
Figure 1-3. 

Pierce (ref. 2) refers to the cusped intersection of the 
two caustic sheets in Figure 1-3 as an "arete."* Beyond the arete, 
which is also called the "focus" or "caustic cusp" in this report, 
the shock folds over upon itself. The fold-over is confined to the 
region bounded by the two caustic sheets. The above description is 
confined to weak shocks. As will be shown in Chapter 5, the fold- 
over occurs for weak shocks, but not for strong shocks which tend 
rather to straighten out smoothly. This verifies a hypothesis of 
Pierce (ref. 16). 

Apparently, there have not been any direct measurements 
of sonic boom focusing caused by local atmospheric inhomogeneities, 
such as a cold spot as illustrated in Figure 1-1, or wind shear fluc- 
tuations as illustrated in Figure 1-4. On the other hand, there have 
been some measurements of focusing generated by maneuvering aircraft 
in supersonic flight. Schematic drawings of shock focusing resulting 
from turning and diving maneuvers are shown in Figures 1-5 and 1-6, 
respectively. Both Figures 1-5 and 1-6 have been redrawn from refer- 
ence 3. Wanner (ref. 4) reports measured focus factors (defined as 
the overpressure at the focus divided by the nominal overpressure) 
up to about 5 for level turns. Maglieri (ref. 5) reports measured 
focus factors up to about 4 for the same type of maneuver. For the 
case of turn-entry, in which an arete similar to the one in Figure 1-3 
is formed. Wanner reports measured focus factors of about 9. 


* Computations of arete locations are included in the Boeing geo- 
metric acoustics program (ref, 15). 
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Most of the pressure-intensification data that is avail- 
able refers to the so-called "sonic cut-off" phenomenon. The cut-off 
phenomenon occurs when an aircraft is flying faster than ambient sound 
speed but slower than sound speed at the ground. At the altitude at 
which the speed of sound is equal to the speed of the aircraft, the 
down-going wave front is reflected into an up-going wave front. The 
locus of points at which reflection occurs represents a caustic sheet. 
Figure 1-7 illustrates the situation for an accelerating aircraft, in 
which the altitude of the point of reflection moves down toward the 
ground with time. 

It should be emphasized that the sonic cut-off intensifica- 
tion phenomenon, which produces overpressure intensification factors of 
the order of 2 (see Maglieri et al . , ref. 6), is to be distinguished 
from proper focusing. In proper focusing, which is the phenomenon of 
interest here, much higher intensification factors are expected. 

Confusing terminology exists in the literature. The cusp 
of an N-wave shock at a single caustic surface has been called a "focus" 
and the associated overpressure a "superboom." When the shock cusp 
occurs where two caustic sheets meet (at the "arete"), the shock cusp 
has been called a "super-focus" and the associated overpressure a 
"super-superboom," respectively. The distinction between the single 
caustic surface and the cusped caustic surface has been discussed by 
A. Pierce (ref. 2). 

Except for Pierce's scaling law analysis (ref. 2), all the 
previous theoretical investigations of sonic boom intensification have 
been confined to phenomena associated with smooth single caustics, i.e., 
superbooms. The analyses of Hayes (ref. 7) and Seebass et al . (ref. 8) 
fall into this category. In contrast, the present report is concerned 
with proper focusing of a shock at a cusped double caustic, i.e., super- 
superbooms. 

The primary approach adopted in this report is different 
from those of the past. Rather than attempt a correction to geometric- 
acoustic theory, or make restrictive assumptions about the nature of 
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the flow field, the full inviscid conservation equations are retained. 
A numerical solution is obtained through the use of a moving mesh that 
propagates with the N-wave (see Chapter 3). By confining the mesh to 
the spatial region of interest, the computer problems which would be 
encountered in following an N-wave for large distances with a code 
using a finite stationary grid are avoided. 

A secondary approach based on Whitham's approximate model 
(ref. 9) of shock wave propagation is also presented. The Whitham- 
type model deals with a single shock wave, whereas the more rigorous 
model described in Chapter 3 treats the focusing of the entire N-wave. 

The contents of the remainder of the report are summarized 

as follows: 

Chapter 3 

The shock following code called GODUNOV is described. 
GODUNOV computes the flow field within the N-wave as it propagates 
through the focus.. The technique employed to solve the conservation 
equations is discussed, as well as the results of some test problems. 

Chapter 4 

The single-shock model code called WHITHAM is described. 
The model is based on a ray-tube-shock-segment formulation in which 
an empirical formula is used to relate ray-tube areas and shock- 
segment Mach numbers. WHITHAM is used as an independent auxiliary 
calculation that can follow the behavior of single curved shocks in 
much less computer time than would be required with GODUNOV. 

Chapters 5 and 6 

Solutions have been obtained for focusing problems involv- 
ing two types of assumptions, namely, 

(a) initially plane N-wave fronts, refracted into concave shapes 
by passing through cold spots 

(b) initial concave-front configurations with prescribed geometric 
parameters such as curvature and rate-of-change of curvature. 

In both cases, the solutions are carried through the focus. Case (a) 
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represents the full problem starting from the physically expected 
initial condition. However, Case (b) is useful for the following 
reasons. First, a significant amount of computer time is saved by 
avoiding the early part of the calculation, namely, that dealing with 
the propagation through the cold spot. Second, having geometric con- 
figurations based on analytic formulas (such as Gaussian and polyno- 
mial functions) allows one to study scaling laws involving, for exampl 
the curvature and the derivatives of the curvature. Third, generally 
prescribed configurations are applicable, not only to cold-spot refrac 
tion, but to the refraction produced by any of several possible physi- 
cal mechanisms. Hence, it is understood that, for the problems involv 
ing assumption (b), atmospheric disturbances or aircraft maneuvers can 
produce the refracted fronts. The flow variables as functions of two- 
dimensional space and time are obtained. Tests are made with changes 
in numerical parameters such as numbers of grid points in order to 
obtain numerical ly-convergent solutions. Focus factors are given for 
various values of the physical parameters. Also investigated is the 
question: Under what conditions will a concave shock fold over? 

Chapter 7 

Some conclusions that follow from our results are summar- 
ized. These include the confirmation of the geometric-acoustics wave- 
folding phenomenon for weak shocks, (Ap/p « 1) as well as the absence 
of wave-folding predicted by Whitham for strong shocks (Ap/p ^ l). 
Focus factors of 19 and 13 have been obtained for polynomial and 
Gaussian front N-waves, respectively. In both cases, high over- 
pressures are confined to spatial regions with scale lengths of the 
order of the wavelength. These results demonstrate the capability 
of solving the focusing problem with the numerical hydrodynamics formu 
lation described in the report. Computations with initial conditions 
representative of proposed supersonic transport operations can now be 
carried out. 
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2. SYMBOLS 


a 

A 


a 

cross 

B 

c 

e 

E 

i 


j 

K 

£ i 

m 

M 

P 

r 


Ar i 


t 

u 


i 


i »J 


= Ar./Ar in GODUNOV 

= area of cell boundary in GODUNOV, area of shock segment 
in WHITHAM 

= cross sectional area of cell in r,z plane 
= normal component of cell boundary velocity 
= sound speed 

= internal energy per unit mass 

1 2 2 

= total energy per unit mass E = e + j (u + v ) 

= index for horizontal rows in GODUNOV mesh; index for shock 
segments and nodes in WHITHAM 

= index for columns parallel to shocks in GODUNOV 

= constant in area versus Mach number relation used in WHITHAM 

= length of shock segment i (GODUNOV and WHITHAM) 

= mass flux across a wave appearing in Riemann problems 

= Mach number 

= pressure 

= shock velocity of segment i (GODUNOV and WHITHAM) 

= radial coordinate 

= vertical component of shock node velocity in WHITHAM 
= separation of horizontal grid lines in GODUNOV 
= ti me 

= velocity component in z (or x) direction 
= shock node velocity in WHITHAM 
= node velocity in GODUNOV (see Figure 3-2) 


9 


u 


= normal component of velocity flowing into a cell in 
GODUNOV 

v = velocity component In r (or y) direction 
V = cell volume in GODUNOV 

V w = wave propagation speed 

x = coordinate parallel to direction of propagation 

y = coordinate normal to propagation direction 

z = axial coordinate 

z.j = axial component of shock node velocity in WHITHAM 
y = ratio of specific heats 

6 = symbol appearing in conservation equations. 

6=0 for x,y geometry, 6=1 for r,z geometry 

6. = angle of inclination of segment i to vertical 

p = density 

north, south, east, west cell boundary 
value at beginning of time step 
value at end of time step 
regions 1, 2, 3, 4 in Figure 3-3 
cold spot 

contact surface (see Figure 3-3) 
normal to cell boundary 
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3. THE GODUNOV CODE 


A computer code, called GODUNOV, has been developed to 
solve the full set of two-dimensional time dependent conservation 
equations for the case of a focusing N-wave. The numerical method 
that is employed in GODUNOV is a modification of a technique originally 
devised by Godunov, et al. (ref. 10) to study the shock layer adjacent 
to a blunt body in supersonic flight. Since Godunov's original pre- 
sentation, the Godunov technique has been applied successfully to a 
variety of blunt body problems, e.g. references 11 and 12. To the 
authors' knowledge, this is the first time the technique has been 
extended to a propagation problem. 

One great advantage of the Godunov scheme is that it pre- 
serves the discontinuity across shock waves of arbitrary strength. 

In this respect, it is superior to the standard finite difference 
codes, such as SHELL, which use artificial viscosity to spread a shock 
over several mesh points and tend to obliterate weak shocks. It should 
be pointed out that GODUNOV also treats any internal discontinuity which 
may arise within the N-wave, through its intrinsic artificial viscosity. 

3.1 Mesh Motion and Geometry 

The mesh geometry employed in GODUNOV is illustrated in 
Figure 3-1. The leading and trailing shocks in the N-wave are shown 
as solid lines. The grid points lie on horizontal lines with fixed 
spacing in the vertical direction. 

Since the high pressure region that results from focusing 
does not extend far from the axis of symmetry, it is desirable to 
place most of the grid points near the axis of symmetry. This can be 
achieved by placing the horizontal grid lines close together near the 
axis of symmetry and further apart at large radial distances. The 
separation between horizontal grid lines in GODUNOV follows the geo- 
metric progression 
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a 


Ar i 


Ar 


i-1 


where a is a constant. We have found that the optimum value of a is 
1.05 for these problems. 

In order to follow the N-wave as it propagates, the grid 
points are allowed to move horizontally but not vertically. The N 
grid points within the N-wave {N = 5 in Figure 3-1) move in such a 
fashion that they are always equally spaced between the leading and 
trailing shocks. There are M grid points behind the N-wave (M = 2 in 
Figure 3-1), and they move so as to be equally spaced between the 
fixed left boundary of the grid and the trailing shock of the N-wave. 
The grid points behind the N-wave provide an indication of the net 
disturbance imparted to the atmosphere after the passage of the N-wave. 

The procedure for moving the grid points is the following. 
First the propagation velocity of each segment of the leading and 
trailing shock is computed by solving a one-dimensional Riemann problem 
as described in Section 3.3 (below). Then the projection of each shock 
segment's normal velocity along the x axis is calculated. The propa- 
gation velocity of a node on the leading or trailing shock is deter- 
mined by using an inverse length weighting of the projected velocity 
of the two adjacent shock segments. Using the notation indicated in 
Figure 3-2, the formula for the shock node velocity is 

1 

u i ,JF0RE {i. + £. +1 ) 

A node falling on the contour labeled j between the lead- 
ing and trailing shocks, indicated by the dashed curve in Figure 3-2, 
is given a velocity 


S'+l 


H+l 


sin e. 


sin e 


i+1 J 


( 1 ) 
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u i,j = M i,JAFT + W i, JFORE " U i,JAFT 5 [JFORE - A JAFT) ^ 

where JFORE and JAFT are the j -Indices of the leading and trailing 
shocks, respectively. 

An attractive feature of this floating mesh scheme is that 
it confines the grid points to the continually changing region of 
interest. Thus, for a given number of grid points, it allows for a 
higher resolution of the flow field than codes with fixed Eulerian or 
Lagrangian meshes. 


3.2 Conservation Equations 

The conservation equations which describe the two- 
dimensional unsteady flow of an inviscid fluid are given below. 


Mass 


1_(P U ) 

+ 9 (pv) _ _ 6 PV 

at ax 

ay y 

x-Momentum 

l_(pu) + 1_(P + 

p 

pu ) + 9 (puv) _ <5 puv 

at ax 

9y y 

y-Momentum 

a_( P v) a_{ P uv) 

+ 9 (P + PV 2 ) _ 6 pV 2 

at ax 

9y y 

Energy 

^ [pe t | (u 2 + v 2 )] + |j u [ 

P + P e + \ (u 2 + v 2 )] 

+ |y v [p + P e + | (u 2 + v 2 )] 

= - 6 j [p + P e + | (u 2 + v 2 )] 


(3) 


(4) 


(5) 


( 6 ) 
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In the above equations 6 = 0 for x, y geometry and 5 = 1 
for r,z geometry. These equations together with the equation of state 
represent a set of five nonlinear equations for the five unknowns 
p, p, e, u, and v. A perfect gas law equation of state has been 
incorporated into GODUNOV, i.e. 

p = (y - 1) P e (7) 


In applying the conservation equations to the moving mesh 
in GODUNOV, Eqs. (2) - (6) are integrated over a cell volume (V = cell 
volume). After applying Green's theorem, the result is 


Mass 


A (p V) 
At 


■If 


P A (U - B) 


( 8 ) 


NSWE 


x -Momentum 


A = (p a sin 0 ) w “ (p a sin 0 )e 


y - Momentum 

A = (pAcos e) F - (pA cos e). 


+ X [ puA (u ■ b) ] (9) 


NSEW 


(P A )t, 


+ (pA) s + If pvA (U - B) 


+ 52ttA (10) 

< cross 


NSEW 
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Energy 


4 (pEV) . 
At 


U 

NSEW 


pAE (U - B) + pUA 


3 


Cn) 


The subscripts in Eqs. (8) - (11) refer to the north, south, 
east, and west boundaries of the cell in question (see Figure 3-1). 

The areas of the cell boundaries are denoted by A, whereas A 

cross 

denotes the cross sectional area of the cell in the r,z plane. U 
represents the normal component of inward-flowing velocity and B the 
normal component of the cell boundary velocity. The angle 0 is measured 
from the positive x axis to the east or west cell boundary as illustra- 
ted in Figure 3-1. 

During the course of a time step, the cell volume changes 
as a result of the mesh motion. The new cell volume at the end of a 
time step, V new , must first be calculated before Eqs. (8) - (11) can 
be utilized. The resulting equations used to update the flow vari- 
able are 


(PV)_ = (pV) Qld + At 


'new 


X.U ( u - B >] 


NSEW 


( 12 ) 


(upV) 


new 


( p uV) old + ^ 

- (pA sin e)^ 



NSEW 
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Equations (12) - (15) are used to solve for p new , u new , v new , and 
E new’ res P ectivel y- 

3.3 Riemann Problems 

Before the right-hand sides of Eqs. (12) - (15) can be 
evaluated, the values of the flow variables at the cell boundaries 
must be determined. This is accomplished in GODUNOV by solving a 
Riemann problem across the appropriate cell boundary. 

The Riemann problem describes how an initial discontinuity 
between two uniform regions evolves with time. In this case, the two 
uniform regions are two adjacent cells separated by a cell boundary. 
There are four Riemann problems associated with each cell. Two of 
them involve moving boundaries (east and west), and two stationary 
boundaries (north and south). 

Consider two adjacent cells in the same horizontal row as 
illustrated at the top of Figure 3-3. The boundary between the cells 
has a normal velocity, B, which is calculated by averaging the normal 
components of the node velocities at both ends of the boundary, i.e. 


n _ si ri e 
B = — 2— 




(16) 
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Once the fluid velocities in the two adjacent cells are resolved into 
components normal and parallel to the cell boundary, the solution of 
the one-dimensional Riemann problem proceeds as indicated in the x-t 
diagram in Figure 3-3. 

In general, a compression or rarefaction wave, wave 1, 
will propagate into region (1), and another compression or rarefac- 
tion wave emanating from the interface, wave 4, will propagate into 
region (4). Between wave 1 and wave 4, a contact surface exists. 

To the extent that the width of the rarefaction waves can be neglected 
(acoustic limit), the three discontinuities divide the x-t plane into 
4 uniform regions, labeled (1) - (4) in Figure 3-3. 

At the fore and aft (or leading and trailing) shocks, 
wave 4 is a shock wave which coincides with the cell boundary in the 
x-t plane. At an interior east or west cell boundary, the boundary 
position has a velocity dx/dt = B as computed above, whereas a north 
or south boundary is fixed and its position coincides with the t axis 
in the x-t plane. In each case, the flow properties at the cell 
boundary are set equal to the flow properties of the region in which 
it lies. 

If wave 4 is a shock wave, the Rankine-Hugoniot relations 
dictate that the mass flux across wave 4 is 

m 4 = ~\f + P3 + (y ~ !) P4] 07) 

and, from momentum conservation, 

P 4 " P 3 + m 4 (u n4 " = 0 (18) 
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where u n is the fluid velocity relative to the wave. In the case 
that wave 1 is a rarefaction wave, the one-dimensional unsteady 
isentropic relations provide that 


u 


nl 



"Y yP 1 /p 1 



(19) 


Since the pressure and normal component of fluid velocity 

are preserved across a contact surface, p 2 = p 3 = p cs and 

u n2 = u n3 = U cs' Ec l uatlons (18) and (19) may now be solved for p c$ 

and U to qive 
cs 3 


m l P 4 + m 4 P l + m I m 4 (u nl " W 


cs 


+ m 4 


( 20 ) 


and 


where 


p l - p 4 + m 4 u p4 + m,u 


1 nl 


cs 


l + m 4 


(21) 


m 


1 


(y - 1 ) 
2y 





(Pcs/Pi) 


/p l 

TFnm 


( 22 ) 


An iterative solution for m 4 using Eqs. (17), (20), and 
(22) is obtained at each fore and aft shock segment. The shock 
velocity is then computed as 
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(23) 


m 4 

’n ' %4 + ^ 

This procedure ts simplified at interior cell boundaries 
by using the weak wave relations across waves 1 and 4. In the 
acoustic limit, 


m 


1 


a nd m 4 



(24) 


and no iterative solution is required for p cs and U cs , 

To determine the properties of an interior cell boundary, 
the boundary velocity, B, is compared to the wave speeds, where 


V * u nl - + (U « ' Unl) I 

V w4 = u n4 + 'V TP 4 /p 4 + ^ 4 P (U cs * U n4 ) 


For example, if < B < U cs> then the cell boundary lies in region 
(2) in the x-t plane, and the flow variables at the boundary are: 


P = P cs 
U = U cs 

p = p 2 = P 1 Cp cs /p l )1/Y 



(26) 
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If the cell boundary falls in region (3), the first two relations 
in (26) would remain the same while 


p = p 3 = p 4 ( P C s / P4 )1/T 


(27) 


3.4 Time Step Computation 

The time step used in updating the flow variables is com- 
puted on the basis of a Courant type stability criterion. This cri- 
terion limits the time step to less than the time required for an 
acoustic signal to travel across any cell, in either the x or y 
directions. 

An acoustic signal travels northward across a cell at 
speed c + v, and southward at a speed c - v. Hence the vertical 
direction time step is 


At y ~ max (c + v, c - v) 

In calculating the horizontal direction time step, the motion of the 
cell boundary must be included. An acoustic signal propagates east- 
ward at a speed c + u, and westward across a cell at speed c - u. A 
signal emanating from the west cell boundary will have to travel a 
distance az . + B.-At, where az . is the shortest horizontal leg of 

the trapezoidal cell, before it encounters the east boundary. Simi- 
larly, a signal leaving the east boundary will travel a distance 

az . - B At across the cell. Therefore, the horizontal time step is 

mi n w 
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( 29 ) 


, j. _ I "" 1 

flt x ” max [c t u - B^, c -' u '+ B^) 

In order to avoid difficulties that could arise when the cell bound- 
aries and/or the fluid are moving much faster than sound speed, 
GODUNOV uses a modified horizontal time step of 


AZ 


At x max (c + |u| + |B, 


min 


c + | u | + |BJ) 


(30) 


Following Godunov, et al. (ref. 10), the overall time step 
for a cell is 


At 


4t x At y 
4t x + 4t y 


(31) 


and the time step used during a cycle is the minimum value of At 
computed for every cell in the grid. 

3.5 Boundary Conditions 

Zero gradient boundary conditions are imposed at the top 
and left boundaries of the grid. This is implemented in GODUNOV by 
setting the flow variables at the west boundary of a cell in the 
first (far west) column of the grid equal to the corresponding values 
at the center of the cell. Similarly, the north boundary of a cell 
in the top (far north) row of the grid is assumed to possess the same 
properties as the center of the cell. 

The right boundary of the grid coincides with the leading 
shock of the N-wave. The appropriate Riemann problem is solved (as 
described in the previous section) for each shock segment in order 
to compute the jump conditions across the shock. The flow field 
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ahead of the shock is assumed to be undisturbed except for the pres- 
ence of the cold spot. 

A symmetry condition exists at the bottom boundary of the 
grid. In the early stages of development, the symmetry was accounted 
for by placing the center of the first row of cells on the x axis, 
and constraining the east and west cell boundaries on the first row 
to be vertical. This configuration is illustrated at the top of 
Figure 3-4. The vertical component of velocity in the first row of 
cells must be zero in this arrangement. 

Several runs with this symmetry configuration produced at 
late times an abrupt change in the slope of the shock between the 
first and second shock segments. In other words, both shocks in the 
N-wave tended to be inclined upon passing through the focus, and the 
constraint of a vertical segment on the axis was artificial. Down- 
stream of the focus the discontinuity in slope seemed to propagate 
upward along the shocks, and it appeared as if the shocks "broke up." 
This situation is illustrated in Figure 3-5, and is a numerical 
artifact. 

To alleviate this "break-up" the symmetry condition was 
re-posed in terms of an imaginary row of cells across the axis of 
symmetry. As the bottom of Figure 3-4 indicates, the imaginary row 
of cells was taken to be the mirror image of the first row. Now the 
first shock segment can be inclined and a vertical component of 
velocity is allowed in the first row of cells. Subsequent runs with 
this configuration produced smooth shock profiles. 

3.6 Test Runs 

A series of test runs was conducted with GODUNOV before 
it was used for the two-dimensional N-wave focusing problem. 

The first test case concerned the one-dimensional propa- 
gation of an N-wave into a uniform atmosphere, without refraction. 
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The Initial relative overpressure of the N-wave was A P 0 /P 0 = 0.1, and 
the initial wavelength was 50 meters. Ten axial grid points were 
placed within the N-wave and two behind it. Figure 3-6 is a plot of 
the pressure distribution both initially and after it has propagated 
for 0.1807 secs (50 time steps or cycles). 

The half wavelength, z, of the N-wave should increase 
with time according to the formula (ref. 13) 



where £ and Av q are the initial half wavelength and velocity ampli- 
tude, respectively. For the parameters of the N-wave in Figure 3-6, 
(av q = 24.34 meters/sec), Eq. (32) predicts a value of z/i Q = 1.100 
after 0.1807 secs. The value computed in GODUNOV was Z/Z Q = 1.090, 
which is within 1 percent of the theoretical value, indicating excel- 
lent agreement. 

The same problem was run on the SHELL code, which is one 
of the standard hydrodynamic codes with a stationary mesh. The pres- 
sure profiles obtained with SHELL are illustrated in Figure 3-7. The 
artificial viscosity in SHELL has spread the shock waves to such an 
extent that they are barely recognizable as discontinuities with well- 
defined amplitudes. It is apparent that SHELL is not capable of 
following shock waves with relative overpressures much less than 0.1, 
which is the range in which we are interested. 

The other test problem run on both SHELL and GODUNOV was 
a numerical simulation of a "cylindrical shock tube" problem. In 
this problem, the ordinary planar diaphragm separating the high and 
low pressure gases is replaced by an imaginary cylindrical diaphragm. 
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At time t = 0, the pressure ratio across the diaphragm is 4.45, the 
diaphragm is instantaneously removed, and the initial discontinuity 
resolves itself into a shock wave and a rarefaction wave propagating 
in opposite directions. The resulting flow is in the radial direction 
only. 

The results of computations with SHELL and GODUNOV are 
illustrated in Figure 3-8. The initial pressure across the diaphragm 
(4.45) was chosen so as to produce a shock wave with a relative over- 
pressure of 1 in the axial flow case. There is no equivalent analyti- 
cal solution to the radial flow problem, but both SHELL and GODUNOV 
indicate that the shock wave is slightly weaker than it would be for 
axial flow. The excellent agreement between SHELL and GODUNOV con- 
firms that GODUNOV is computing the radial flow correctly. This con- 
firmation together with the axial flow test case results verifies that 
GODUNOV is a valid two-space-dimensional fluid dynamic code. The fact 
that the shock discontinuity is smeared slightly more by GODUNOV than 
by SHELL (in Figure 3-8) shows that the intrinsic artificial viscosity 
is slightly greater in GODUNOV than in SHELL. 
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4. THE WHITHAM CODE 

An auxiliary computer code, entitled WHITHAM, has been 
created as a supplement to GODUNOV. WHITHAM follows the propagation 
of a single shock according to the approach described by Whitham in 
reference 9. In our case, the shock represents the leading shock of 
the N-wave, and the validity of the WHITHAM code is contingent upon 
a lack of interaction between the shock and the flow behind it. 

Whitham' s original premise was that a curved shock may be 
envisioned as a chain of planar elements, each of which propagates 
down a tube of varying cross section. The propagation velocity is 
determined by an empirical area versus Mach number relation. Thus, 
as the shock begins to focus, the segments near the focus are "com- 
pressed" and their propagation velocity increases. In this respect, 
WHITHAM is a higher order formulation than the ordinary geometric 
acoustic ray tube concept where every point on the wave front propa- 
gates at the local sound speed. 

4.1 Numerical Description 

A sketch of the shock front in WHITHAM is illustrated in 
Figure 4-1. Each line segment, which represents in the figure the 
cross-section of a planar element, has a velocity, q.,- normal to 
itself, where is equal to the segment Mach number, M.. , multiplied 
by the ambient sound speed. The differential relationship between 
Mach number (M) and "segment area" (A) is the one proposed by Whitham 
(ref. 9), i.e. 


dA = - 2MdM 

A (M 2 - 1) K(M) 


(33) 


where K(M) is a slowly varying function of Mach number, given in 
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reference 9. Since K.(M) only ranges from 0.5 for weak shocks to 0.394 
for strong shocks, Wh.tth.am (ref. 9) suggests that Its variation may be 
neglected In Integrating Eq. (33), The resulting simple relationship 


A^ (M^ - 1) = const (34) 

is the one employed in WHITHAM. The area A^ associated with each 
shock segment is (a) the area of revolution of the segment about the 
axis of symmetry for r,z geometry, or (b) the area of a strip, per 
unit length in the z-direction, in x,y geometry. 

The sequence of events occurring during one cycle, or 
time step, in WHITHAM is the following. 

1. The lengths and areas of the shock segments are computed from 
the r,z (or x,y) coordinates of the shock nodes. 

2. Mach numbers for every shock segment are computed from Eq. (34) 

(K is read in as input data; we have used K = .5 for most runs), 
yielding segment velocities, q^ . 

3. The time step. At, is computed by taking the smallest value of 
£.j/q.j. This criterion prevents the segments from moving a 
distance larger than their own length in one time step. 

4. Using the shock segment velocities and geometry, compute the 
r and z (or x and y) components of the shock node velocities. 

The method used is described in detail in the next section. 

5. Move each shock node a distance r^At in the radial direction 
and z^At in the axial direction. 

4.2 Shock Node Velocities 

Two different methods have been used to compute the shock 
node velocities, and z^ , from the segment velocities and geometry. 

In the first method, the node velocity components are 


26 



computed as an inverse length weighting of the adjacent segment 
velocities. The equations are 



In the second method, the two adjacent shock segments are 
displaced a distance q^t normal to themselves and the new node posi- 
tion is computed as the geometric intersection of the two segments. 
The equations are 



Equations (37) and (38) are singular when the two segments 
are parallel. Therefore, Eqs. (35) and (36) are used only when the 
cosine of the angle between the segments differs from 1 by more than 
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5. RESULTS FOR INITI ALLY-CONCAVE N-WAVES 


5.1 Polynomial Front 

Two different initially-refracted configurations are 
investigated, using GODUNOV. The first configuration to be discussed 
is called the "polynomial front" N-wave. Both the fore and aft shocks 
of the polynomial -front N-wave at t = 0 are represented by the formula 



0 < r < 100 (39) 


where z and r are in meters, and z Q locates the fore shock and the aft 
shock. 

Consider the 4 pairs of approximately equally-spaced curves 
shown in Figure 5-1. For each pair, the curve on the right represents 
the fore shock profile at a given time, while the curve on the left 
represents the aft shock profile at the same time. The first pair of 
curves (labeled t = 0) corresponds to a polynomial -front N-wave with 
an initial wavelength of 10 meters and a radius of curvature at the 
axis of 50 meters. The initial relative overpressure is chosen to be 

_3 

(Ap/p) o = 10 , where Ap is the pressure jump across the fore shock 

and p Q is the ambient pressure. 

The shock profiles at the four indicated times in Figure 
5-1 illustrate how the N-wave changes shape as it propagates to the 
right into a uniform atmosphere. Note that the point of inflection 
on either front migrates toward the axis as the N-wave approaches the 
focus, defined here to be the position of maximum pressure, which is 
located at z = 79 meters (corresponding to t = .172 sec, not shown in 
Figure 5-1). According to geometric acoustics, the focus is located 
at the center of curvature of the fore shock, at z = 70 meters in 
this case. We will designate the geometric-acoustic focal point as 
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the "nominal focus," to distinguish it from the "actual focus" deduced 
from the calculation. At the actual focus, the shock, develops a cusp 
at the axis of symmetry. The profile at the focus resembles the geo- 
metric-acoustic wave-folding picture in Figure 1-3. 

The incipient break-up of the last shock profiles at the 
axis in Figure 5-1 is a numerical artifact involving use of a verti- 
cal shock segment at the axis. A suitable modification which avoids 
this break-up has been discussed in Section 3.5. 

The polynomial front N-wave shown in Figure 5-1 and 
defined above was run several times with different grid sizes. 

Although the shock profiles exhibit the general shapes shown in 
Figure 5-1 for all of the grids, the focus factors, A P max / A P 0 > depend 
on the number of grid points used. Table 5-1 displays the focus 
factors corresponding to each of the grids. The authors believe that 
the grid consisting of 100 equally spaced radial points and 20 axial 
points yields a reasonably accurate value for the focus factor, namely 
18.7, in the sense that increasing further the number of grid points 
will not change this value significantly. 

Moreover, the use of inclined cell segments on the axis 
is believed to yield more accurate results (see Section 3.5). Through- 
out the rest of this chapter, the results refer to inclined cell seg- 
ments. Focus factors obtained with inclined cell segments are 
approximately 10 percent higher than those obtained with vertical 
cell segments. Thus, a better estimate of the focus factor for the 
polynomial front is about 20. It is interesting to note that the 
focus factor is 11.5 at the position of the "nominal" focus (geo- 
metric acoustics). 

Figure 5-2 illustrates the pressure profiles (as computed 
with the 100 x 20 grid) along the axis at three different times. The 
last curve in Figure 5-2 corresponds to the pressure signature at the 
focus. The formation of spikes near the- fore and aft shocks of 
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the N-wave Is evident. Detailed structure behind the aft shock, is not 
shown because the zoning was coarse behind the N-wave. 

The pressure profiles at three different radial positions 
are illustrated in Figure 5-3. The lowest curve in Figure 5-3, which 
is the profile on the axis, is the same as the last curve in Figure 
5-2, drawn to a different scale. The rapid drop of pressure with 
distance from the axis is apparent in Figure 5-3. 

Figure 5-4 is a plot of the relative overpressure behind 
the fore shock of the N-wave versus axial position of the shock. 

The solid curve represents the results of the 100 x 20 grid GODUNOV 
calculation, while the dashed curve represents the equivalent (100- 
point) WHITHAM calculation. The maximum relative overpressure com- 
puted with WHITHAM (.0285) is 52% higher than the maximum relative 
overpressure computed with GODUNOV (.0187). 

This is apparently due to the fact that there is no 
rarefaction wave in WHITHAM to relieve the pressure buildup. The 
calculations in WHITHAM are terminated when the shock segments 
overlap, or cross each other. (This occurs somewhat earlier than 
the GODUNOV focus.) It can be seen from Figure 5-4 that, along the 
axis, relative overpressures greater than .002 (twice nominal) 
occur in a spatial interval 40 meters long; and relative overpressures 
greater than .01 (about half the maximum) occur in a spatial interval 
10 meters long. In the radial direction, the corresponding interval 
lengths are 20 and 10 meters, respectively. 

In order to assess the relative intensities of point foci 
and line foci for the same initial conditions, the problem defined 
above was recomputed with WHITHAM in a two-dimensional x,y geometry 
(line focus). A maximum relative overpressure of .0117 was computed, 
as compared to the value of .0285 for the point focus. 
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The effect of varying the initial strength of the poly- 
nomial shock front was also investigated with WHITHAM . Calculations 
with initial relative overpressures , Ap o /p Q , of 0.01 and 0.10, in 
r,z geometry, resulted in maximum relative overpressures of 0.154 and 
0.63, respectively. These results indicate that the intensification 
due to focusing of stronger shocks is not as severe as for weaker 
shocks. Furthermore, the geometric-acoustics wave-folding description 
of focusing is not valid for shocks with relative overpressures of 
the order of 1 or higher. This is discussed below in connection with 
some GODUNOV calculations for "Gaussian front" N-waves. 

5.2 Gaussian Front 

We define a "Gaussian front" shock wave as one with a 
profile satisfying the equation 



where r and z are both in meters. The profile labeled "t = 0" in 
Figure 5-5 is an example of a Gaussian front. (Only the fore shock 
profiles are shown.) The other curves in Figure 5-5 represent the 
fore shock profile of the N-wave at later times. Initially, the 
N-wave has a nominal relative overpressure Ap Q /p 0 = 10 and a 
wavelength of 10 meters. The formation of a cusp at the axis as the 
focus is approached in Figure 5-5, and the gradual decay of the cusped 
portion of the shock beyond the focus, confirms the geometric-acoustics 
description of the primary shock shown in Figure 1-3. 

The geometric-acoustics wave-folding picture in Figure 1-3 
indicates that the primary shock should be reflected from the axis of 
symmetry. Further, the reflected shock ends at the caustic sheet, and 
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a third (logarithmic) discontinuity joins the "ends" of the reflected 
shock. The question arises whether this structure may be inferred 
from our calculations. The computed two-dimensional pressure distri- 
bution immediately behind the fore shock just after it has passed the 
focus has been plotted in Figure 5-6. A portion of the mesh used in 
the computation is shown in Figure 5-6, and the numbers within the 
cells represent the pressure at that location. If the presence of a 
reflected shock and a "logarithmic discontinuity" are to be inferred 
fromthe computations, there should be a jump in pressure as one scans 
from right to left along a horizontal row in the mesh. This type of 
pressure jump does not appear in Figure 5-6. It is perhaps not sur- 
prising that the GODUNOV calculations do not reveal the presence of 
secondary discontinuities, because the strong rarefaction behind the 
fore shock probably swamps such discontinuities. As opposed to this, 
the wave-folding picture shown is for a single shock without a strong 
rarefaction behind it. 

A tabulation of focus factors calculated with different 
numbers of mesh points for the Gaussian-front N-wave discussed above 
is shown in Table 5-2. The relatively small change in focus factors 
between the last two grids in Table 5-2 indicates that approximate 
convergence has been obtained. The results discussed in this section 
refer to the 50 x 50 grid, which required 34 minutes of CDC 6600 
computer time. 

Pressure signatures along the axis of symmetry at three 
different times are plotted in Figure 5-7. Here again, the develop- 
ment of steep spikes adjacent to the fore and aft shocks of the N- 
wave is apparent. The pressure signature at the time corresponding 
to focusing in Figure 5-7 indicates a focus factor of 13.0 (relative 
overpressure = .013) for the Gaussian front N-wave. 

Note that the location of the actual focus (z = 64 meters) 
is close to the center of curvature (the nominal focus) of the initial 
N-wave (z = 70 meters). 
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For the Gaussian-front N-wave, relative overpressures 
larger than .002 (twice nominal) occur in an axial interval 35 
meters long, and relative overpressures greater than .0065 (about 
half the maximum) occur in an axial interval 20 meters long. Tn the 
radial direction, the corresponding interval lengths are 6 meters 
and 1.5 meters, respectively. 

5.3 Strong Shocks 

The same Gaussian front N-wave described above was rerun 
on GODUNOV with a much higher initial overpressure in order to inves- 
tigate the focusing of relatively strong shocks. The shock profiles 
for an N-wave with Ap Q /p o =0.90 are shown in Figure 5-8. Notice the 
lack of a cusp in the shock profile at the focus in this case. In 
fact, the entire picture looks more like Whitham's picture of a con- 
cave shock "overshooting" than the geometric-acoustics wave-folding 
picture (compare figures 1 and 4 in ref. 9). The focus factor is 1.5. 

GODUNOV has been used to investigate the propagation of 
other relatively strong concave and convex shocks. The convex shock 
that was studied is identical to one of the shocks that Collins and 
Chen (ref. 14) used in their study of shock wave diffraction. 

As indicated in Figure 5-9, the initial shock profile is 
composed of three straight sections labeled A, B, C. Segments A and 
C have Mach numbers of 2.23 (Ap Q /p o = 4.63), while the inclined seg- 
ment, B, has a Mach number of 1.576 (Ap Q /p o = 1.73). The decay of 
the convex portion of the shock at later times as illustrated in 
Figure 5-9 is in close agreement to Collins' and Chen's results. 

A similar shock, with a concavity instead of a convexity, 
is shown in Figure 5-10. Although the shock does straighten out, the 
return to a planar configuration is not as smooth as it is for the 
convex shock. The spike in the last shock profile is a numerical 
artifact due to the symmetry constraint employed in that particular 
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run. (The modifications referred to in Section 3-5 of Chapter 3 
would eliminate this artifact.) 

Comparisons of Figure 5-8 and 5-10 with Figures 5-1 and 
5-5 reveal that the focusing of concave shocks is much different for 
weak shocks than it is for strong shocks. Weak shocks with relative 
overpressures much smaller than unity focus according to the geometric- 
acoustic wave-folding mechanism, whereas strong shocks with relative 
overpressures of the order of unity or higher tend to straighten out 
as Whitham (ref. 9) predicted. 
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6. RESULTS FOR COLD-SPOT REFRACTION OF INITIALLY-PLANAR N-WAVES 

In early runs, spherical cold-spots with uniform tempera- 
tures below ambient were investigated. The incident shocks were 
assumed to be initially planar, which is a good approximation if the 
radius of the cold-spot is much smaller than the radius of the Mach 
cone associated with the bow shock of a supersonic aircraft. The 
discontinuity in temperature at the cold-spot interface caused a 
reflected shock. The disturbance in the N-wave (as calculated by 
GODUNOV) that was caused by the reflected shock produced complicated 
solutions without providing further insight into the focusing mechanism. 
Consequently, in later runs a continuous transition in temperature was 
imposed at the cold-spot boundary, so that reflected shocks, if any, 
were weak and did not appear. The temperature variation within the 
cold-spot was taken to be 


T_ 

T 


1 + — exp 
p o 


-3 (z 


z + r^ /r^ 

spot' 'spot 


exp (-3)| (41) 


where z spot is the location of the center of the cold-spot, r $ p Qt is 
the cold-spot radius, and Ap/p Q is the relative density change between 
the center of the cold-spot and ambient conditions. The pressure in 
the cold-spot is taken to be the same as the ambient pressure. 

The results of a run in which the temperature transition 
is continuous are sliown in Figure 6-1. The solid curves represent 
the fore shock and the dashed curves the aft shock. This figure shows 
a numerical break-up occurring at late times, which is an artifact and 
was corrected in later runs (see Figure 6-2). The break-up occurs 
after the focus and has a negligible effect on the value of the focus 
factor. Therefore it may be ignored in the following remarks. 
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In the results illustrated in Figure 6-1, Ap/p Q = 1.0, 

z . = 120 meters, and r . = 50 meters. The first pair of curves 
spot spot r 

represent the impinging planar 'N-wave. The second pair of curves in 
Figure 6-1 illustrate the refraction that is caused by the cold-spot 
slowing down the inner portion of the N-wave. Note the resemblance 
between the cold-spot refracted shock fronts in Figure 6-1 and the 
Gaussian front of Figure 5-5. The initially-refracted concave shocks 
described in Chapter 5 started with uniform overpressures along the 
shock fronts. At the time the N-wave shown in Figure 6-1 emerges from 
the cold-spot, the overpressure variation along the fore shock is 55%, 
i.e., relatively small compared with 390% at the focus (see below). 

The third pair of curves in Figure 6-1 have been drawn at 
a time when the fore shock has already propagated past the focus, which 
occurs at z = 231 meters. With a grid composed of 20 radial points and 
7 axial points, a focus factor of 3.9 was computed for the problem 
shown in Figure 6-1. A finer mesh would' result in a larger focus 
factor; if the results shown in Table 5-1 can be used as a guideline 
to extrapolate to a finer mesh, a focus factor of 11.3 can be estimated. 

GODUNOV has also been employed to compute the two- 
dimensional (x,y geometry) cold-spot focusing that results from a 
situation equivalent to the one shown in Figure 6-1. In other words, 
the cold-spot is now cylindrical instead of spherical so that a line 
focus will result instead of a point focus. A focus factor of 1.5 
was calculated using the same 20 x 7 mesh. Extrapolation to a finer 
mesh in this case would lead to a focus factor of 4.4. 

Figure 6-2 shows the results of a computation with a mesh 
consisting of 50 radial points and 50 axial points. The cold-spot 
parameters are Ap/p q =1.0, z S p 0t = 250 meters, and r spot = 150 meters. 
Thus, the cold-spot in Figure 6-2 is larger than the one in Figure 6-1, 
although the central temperatures are the same. The three pairs of 
curves in Figure-6-2 represent the N-wave (a) impinging on the cold- 
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spot, (b) midway through the cold-spot, and (c) just after focusing 
(the focus is located at z = 430 meters). Note the absence of a 
break-up in the last pair of shocks in Figure 6-2. This is due to 
the improved symmetry condition employed in this calculation (see 
Section 3.5). 

A plot of the relative overpressure behind the fore shock 
versus axial position of the fore shock is given in Figure 6-3. Sig- 
nificant increases in the overpressure are not observed until the 
fore shock is almost through the cold-spot. This is due to the small 
time lag between refraction and focusing. The focus factor in this 
calculation is 16.7, and the axial distance from the focus at which 
the overpressure is one-half the maximum value is about 30 meters 
(.6 wavelengths in this case). 

All of the results quoted above are for cold-spots with 
center temperatures equal to one-half the ambient temperature. This 
is an example of an extreme temperature inhomogeneity that would not 
be encountered in the real atmosphere. Less extreme temperature 
inhomogeneities would produce lower focus factors. However, other 
inhomogeneities, such as wind shear fluctuations, may be more signi- 
ficant but were not considered in the calculations. 
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7. CONCLUSIONS 


The results of this investigation demonstrate the capa- 
bility of GODUNOV as a two-space-dimensional shock-following code for 
calculating sonic boom N-wave focusing. 

The flow field at and near the focus has been computed 
by GODUNOV for two types of sample problems. The first type of prob- 
lem involves the refraction and subsequent focusing of a planar N-wave 
by a cold-spot. Focus factors of 11.3 and 16.7 are obtained for the 
two spherical cold-spots (point foci) investigated. A cylindrical 
cold-spot (line focus) similar to the first spherical cold-spot yields 
a focus factor of 4.4. 

The second type of focusing problem investigated in this 
study concerns the focusing of an N-wave with a concave front of pre- 
scribed shape. The curved front might be caused, for example, by 
atmospheric refraction or by aircraft maneuvers. Two different initial 
shock-front shapes are studied. In one case we obtain a focus factor 
of 13, and in the other case a factor of 20. In all cases, the focal 
region (as defined by the distance from the focus at which the N-wave 
overpressure falls to one-half the maximum overpressure) extends no 
more than 3 wavelengths from the focus. These results illustrate the 
dependence of the focus factor on the initial shape of the shock front. 

The study has also provided valuable insight into the 
process of focusing. The wave-folding mechanism predicted by geo- 
metric acoustics (Figure 1-3) for a concave shock has been verified 
for weak shocks, although no evidence of secondary discontinuities, 
i.e. reflected shocks, has been observed. Wave-folding is the mechanism 
responsible for sonic boom focusing. On the other hand, strong concave 
shocks with relative overpressures , Ap o /p Q , of the order of unity or 
higher tend to straighten out or overshoot rather than fold over. 

Thus, the hypothesis of Whitham (ref. 9) that a shock will straighten out 
without fold-over, and the fold-over hypothesis of geometric acoustics 
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(refs. 2, 16), are complementary to one another; both are valid but 
under different pressure conditions. 

It may be noted that the GODUNOV code can be linked to 
a geometric-acoustics code (such as the codes described in references 
1 and 15) in order to calculate sonic boom signatures from maneuvering 
aircraft. Predictions of the intensity and extent of "super-super- 
booms" (see Introduction) from prescribed aircraft maneuvers could 
provide the basis for defining acceptable flight operations for 
supersonic aircraft. 
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APPENDIX A 


USER'S MANUAL FOR PROGRAM GODUNOV 

GODUNOV is designed to compute the flow field that results 
from the focusing of an N-wave of arbitrary strength. The focusing 
is associated with a concavity in the wave front, which can be caused, 
for example, by a cold-spot in the path of the N-wave. One may also 
prescribe the initial concavity in the wave front. Thus, GODUNOV can 
be used to study focusing of (i) an initially straight-front N-wave 
upon passing through a cold-spot of prescribed size and intensity, or 
(ii) a curved-front N-wave of prescribed shape in a homogeneous 
atmosphere. 

The input to GODUNOV determines the initial state of the 
fluid, the mesh spacing, and the frequency of printed output. A 
short output giving key data such as the shock positions and the 
positions and value of maximum pressure in the N-wave, is printed at 
every cycle. A long output describing the full two-dimensional state 
of the fluid is printed at the desired cycle frequency. In addition, 
tape dumps are made periodically to store information for a possible 
restart at a later date. 


Input Data 

The input data required to start a problem consists of a 
six-card package as described below. Any self consistent system of 
units can be employed. The cgs system is used for the problems des- 
cribed in this report, i.e. lengths are in cm, densities are in 

2 

gm/cc, and pressures are in dynes/cm . 


40 



1st CARD 


Columns 1 - 10 


Columns 11 - 20 


Columns 21 - 80 


The problem number (CPRQBl ts required 
in a F10.3 format 

The cycle Ctime step) number CCYCLE) is 
required in a F10.3 format. CYCLE = 0 
to start a new problem. In restarting 
in order to continue a problem, CYCLE 
is the last cycle number of the previous 
run. 

are left blank 


The data in cards 2-6 are input in the Namelist format. 
The Namelist feature provides considerable flexibility by requiring 
only input that specifies the user's choice of options or is different 
from the preset data. However, it is only available on certain com- 
puters, e.g. CDC 6600. The procedure for inputting data via a Name- 
list format can be found in most FORTRAN IV manuals. 


2nd CARD - $PRELIM 

The following variables are contained in Namelist PRELIM: 

IMAX the number of grid points, or rows, in the 

r (radial or vertical) direction, indexed by i. 

JF0RE the number of grid points, or columns, in the 

z (axial or horizontal) direction. JF0RE is 
also the value of the j index corresponding to 
the fore shock. 

JAFT the number of grid points, or columns, in the 

z direction behind and including the aft shock. 
JAFT is also the value of the j index corre- 
sponding to the aft shock. 

TMAX time in seconds at which computations are to 

be terminated. 

CYMAX cycle number at which computations are to be 

terminated. (The program will stop computing 
whenever T > TMAX or CYCLE > CYMAX, whichever 
occurs first.) 
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CPRINT 

cycle interval between long outputs 

CDUMP 

cycle interval between tape dumps 

T 

value of T (time) at start of run 

GAMMA 

ratio of specific heats 

CARD - $MESH1 


Namelist MESH 1 contains the variables: 

RMAX 

value of r at outer boundary of grid 

DZFORE 

initial (uniform) spacing of grid contour 
lines within N-wave and parallel to shocks 

DZAFT 

initial spacing of grid contour lines behind 
N-wave, i.e. between the left boundary and 
the aft shock, and parallel to shocks 

ZSPOT 

value of z coordinate corresponding to center 
of cold-spot 

DZSPOT 

diameter of cross-section of cold-spot 

A 

ratio of spacing between successive radial 
grid lines 

ZAXIS 

array representing z coordinate of grid points 
on axis of symmetry (computed internally for 
standard problem corresponding to OPT = 0, 1, 
or 2 (see below) ) 

OPT 

integer indicating initial shock front shape 
OPT = 0 planar front 
OPT = 1 polynomial front 

z = z Q + RMAX [(r/ RMAX) 2 - \ (r/RMAX) 4 ] 

OPT = 2 Gaussian front 


z = z + R ^} X jjl - exp [-20 (r/RMAX) 2 ] J 


where z locates the fore and aft shocks 
and the contours in between 
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4th CARD - $MESH2 

Namelist MESH2 contains the following variables: 

R(i) vector consisting of tfie r-values of the fixed 

horizontal grid lines 

Z(i,j) two-dimensional array consisting of the z-values 

of the mesh points. The first index (i) labels the 
horizontal line, and the second index (j) labels the 
node on the horizontal line between the shocks. 

When one of the standard problems corresponding to 
OPT = 0, 1, 2 is to be run, MESH2 should be left blank, i.e. 

$MESH2 $ 


5th CARD - $STATE1 

STATE1 allows the initial fluid state to be described in 
compact form. The following variables are in STATE1. 

PINIT(k) 6-component vector specifying fluid state ahead 

of N-wave (see 6th card) 

PINIT(l) pressure 

PINIT(2) density 

(The other 4 components are not used. Energy and sound 
speed are computed in the program, assuming a static state.) 

PSP0T(2) 1 + ratio of density at center of cold-spot to 

ambient density. 


(The other 5 components are not used.) 

The variation of density within the cold-spot obeys the 

equation 


p 


P 0 jl + PSPOT (2 ) • exp - 3 ju - ZSP0T)‘ 
- PSP0T(2) • exp (-3)J 
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where RSPOT is the radius of the cold spot, and ZSPOT is the axial 
position of the center of the cold spot. 

S the initial relative overpressure of the N wave, 

t.e. S = ad /p where ap is the pressure jump 

O O 0 

across the fore shock, and p is the ambient 

K o 

pressure ahead of the N wave. 

The program sets up a linear variation of all the state 
variables within the N-wave. The isentropic relations are used to 
relate the pressure to the other flow variables. The distribution of 
all the flow variables are approximately symmetric about ambient 
conditions (i.e., about the mid-point of the N-wave). 

6th CARD - $STATE2 

Namelist STATE2 allows the user to specify x,y or r,z 
geometry, as well as an arbitrary initial distribution of flow variables. 
The following variables are in STATE2 

IXY IXY > 0 implies x,y geometry 

IXY < 0 implies r,z geometry 

PINIT(k) The ambient flow field ahead of the N-wave can 

be specified as a vector as follows: 

PINIT(l) = ambient pressure 

PINIT (2) = ambient density 

PINIT (3 ) = z-component of ambient velocity 

PINIT(4) = r-component of ambient velocity 

PINIT(5) = ambient total energy per unit volume 

PIN IT (6) = ambient sound speed 


44 



P(i.j,k) 


three-dimensional array describing the initial 
flow field in the N-wave. The first index 
labels the horizontal row of the cell in the 
mesh, the second index (j) labels the cell on 
the horizontal line between the shocks, and 
the third index (k) labels the state variable 
as follows: 

k = 1 - pressure 

2 - density 

3 - z-component of velocity 

4 - r-component of velocity 

5 - energy per unit volume 

6 - sound speed 

To run one of the standard problems, only IXY need be 
input in STATE2. 

The input data package for restarting a problem consists 
of the two cards described below. 


1st CARD 

Columns 1-10 the problem number (CPROB) in a F10.3 

format. 

Columns 11 - 20 the cycle number at which the problem 

is to be restarted in a F10.3 format 


2nd CARD 

Columns 1-10 TMAX, the time at which calculations 

are to be terminated 

Columns 11 - 20 CYMAX, cycle number at which calcula- 

tions are to be terminated 

Columns 21 - 30 CPRINT, cycle frequency for long outputs 

Columns 3.1 - 40 CDUMP, cycle frequency for tape dumps 
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A subroutine-by-suhroutine flow chart of GODUNOV is shown 
in Figure A-l» The purposes of the various subroutines, are as follows. 

INPUT reads all required input data regarding mesh. geometry, 
initial values of the flow variables at all points in space, and physi- 
cal parameters such as the ratio of specific heats. 

NODES deals with changes of the mesh geometry as a func- 
tion of time, and computes (at different points in the logical circuit): 

(a) the lengths of shock segments (connecting adjacent shock nodes) 

(b) the velocities of the shock segments normal to themselves, and 
the corresponding axial velocities of the shock nodes 

(c) the velocities of the moving (non-shock) boundaries of the 
interior cells 

(d) the new positions of all mesh nodes. (The various portions (a, 
b, c, or d) of NODES are called at different points in the flow 
sequence, as indicated in Figure A-l.) 

SHOCK solves the Riemann problem according to the scheme 
suggested by Godunov et al . (ref. 10) at all shock segments (both fore and 
aft shocks), yielding the segment velocities normal to themselves, and 
the values of the flow-variables on both sides of each shock segment. 

These quantities are employed in evaluating the fluxes (of mass, momen- 
tum, and energy) on those special cell boundaries which coincide with 
shock segments (for cells adjacent to the shocks). SHIO ("shock- 
input-output") sets up the input to SHOCK, and processes its output 
for use by the main program. 

RIEMANN solves the Riemann problem at all boundaries of 
each interior cell (not adjacent to shocks), taking into account the 
motion of the moving cell boundaries, yielding the values of the flow 
variables (continuous) on the cell boundaries. RIO ("Riemann-input- 
output") sets up the input to RIEMANN for each cell boundary, and 
processes its output. Accuracy can be maintained while employing a 
linearized version of the Riemann problem for the continuous portions 
of the flow, with a resulting economy in computation time. 
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FLUXES computes the fluxes of mass, momentum, and energy 
at the cell boundaries in preparation to updating the values of the 
flow variables at the cell centers. (FLUXES is called by NEWFLO.) 

NEWFLO computes for a given cell: 

(a) the cell boundary surface areas (of revolution in r,z geometry) 

(b) the old cell volume (of revolution in r,z geometry) prior to 
updating the mesh node positions, and the new cell volume after 
updating the mesh node positions. 

NEWFLO then updates the values of the flow variables at 
the cell center. This latter updating is done by means of the conser- 
vation of mass, momentum, and energy over the volume of the cell. 
First, the currents of mass, momentum, and energy at the cell boundary 
surfaces are obtained by multiplying the fluxes (obtained from FLUXES) 
by the cell boundary surface areas. Summing these currents and multi- 
plying by the time increment gives the change in a quantity Q (not 
indicated in the figure), where Q represents the total mass, momentum, 
or energy contained in the cell. Q is approximated by multiplying 
the volume of the cell by the mass density, momentum density, or 
energy density within the cell (assumed constant across the cell). 

The new value of Q is proportional to the new cell volume. Hence, 
the new value of mass density, momentum density, or energy density 
within the cell is obtained by dividing the new value of Q by the new 
cell volume. The new values of the flow variables within the cell 
are subsequently readily calculated. 

EOS ("equation of state") computes the pressure and sound 
speed when the mass density, energy density, and fluid velocity com- 
ponents are given. We are presently using an ideal gamma-law-gas 
equation. 

DT determines the time interval to be used in updating 
the mesh and the flow variables. The time interval is chosen to sat- 
isfy the Courant stability criterion, which requires that the time 
interval be less than the time required for a sound signal (speed of 
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sound combined with fluid velocity) to travel across any cell, in 
either the radial or axial directions. The motion of the moving 
boundaries must also be taken into account here. 

OUTPUT prints out appropriate information at desired 
intervals of time. 

This completes the description of the subroutines used 

in GODUNOV. 

A listing of GODUNOV follows. 
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RUN VERSION 2.3 


PSR LEVEL 298- 


000003 


PROGRAM GOO UNO V ( INP LT , OUTPUT ,TAPE3»TAPE5=INPUT ,TA PE 6= OUTPUT » T APE 7) 
COMMON P(50,20,6), R(55), Z(55,20), ZSTART (50) * ZEND(50), PINIT(6I 

000003 


1, PSPOT (6) , ZOLO(50> ,PFORE (50,6) ,PA FT 150*69 »FORSEG(50) ,FORVEL<50) , 
2AFTSEG (501 , AFTVEL (5 0) , PN( 20 ,6) , P£(6) , AN(20> ,UFORE(50) , UAFT (50 ) , 
3 VNE (20) * I XT 
4,ZAXIS(20) 

COMMON RN, RS, ZNE, ZNW, ZSE, ZSW, ESEG, HSEG, GAMMA, T, TMAX , OT, 

000003 


1 CYMAX, CPRINT, OTMIN, IDT, JDT , PROB, CYCLE, CDUMP, LREAO, MWRITE 
2, LTAPE, DELR, ISPOT, VOLOLO, VOLNEW, AS, AW, AE , VNORME , VNOPMW , 

3 UNORME, UNORMW, UTANGE, UTANGW, VFORE , V AFT, VNW, VSE,VSW, IMAX , 

4 JFORE , JAFT, I, J, PMAX, HMACHF, HMACHA , IPMAX, JPMAX 

COMMON RSPOT, ZSPOT 

000003 


LREAO = 5 

000004 


MWRITE=6 

000005 


LT APE = 3 

000006 


REWIND LTAPE 

000010 


CALL SETPLTS 

000011 


CALL INPUT 

000012 


IF (CYCLE .GT. 0.) GO TO 5 


C 

START NEW CYCLE 

000015 

1 

IF (CYCLE .GT. 0.) FMAX-0. 

000020 


PMIN = 1.E19 

000022 


00 3 1=1, IMAX 

000023 


CALL NODES 

000024 


CALL SHIO 

000025 


CALL NCDES2 

000026 


DO 2 J=l, JFORE 

000030 


CALL NODES 3 

000031 


CALL RIO 

000032 


CALL OTCALC 

000033 


IF (CYCLE . EQ. 0.) GO TO 2 

000034 


CALL NEWFLO 

000035 


CALL N0DES4 

000036 


CALL NUFLOW 

0000 37 


CALL EOS (GAMMA ,P( I,J,1),P(I,J,2),P(I,J,3),P(I,J,4),P(I,J,5),P(I,J, 

000055 


16) ) 

IF (P(I,J,1) .LT. PMAX) GO Tn 4 

000063 


PMAX = P(I, J,1 > 

000066 


IPMAX = I 

000066 


JPMAX = J 

000067 


4 IF (P(I,J,1) .GT. PMIN) GO TO 2 

000076 


PMIN = P(I,J,1) 

000100 


IPMIN = I 

000101 


JPMIN = J 

000102 


2 CONTINUE 

000105 


3 CONTINUE 

000107 


IF (PMIN .IE. 0.) C YM A X=C YCL F 

000113 


IF (CYCLE , FO . 0.) CTMINO=0. 

000115 


IF (OTMIN .LT. DTMIKC) CT=OTM IN»' 2 70TMIN0 

000121 


IF (OTMIN . 6E . OTMTNC) DT= DT M IN 

000124 


OTMING=OTMIN 

000125 

5 

CALL OUTPUT 

000126 


IF ( AMOD (CYCLE, CPF: TM> .FG. 0. .OP. T . GE . TMAX .OP. CYCLE .GE, 



1 CYMAX) CALL OUT L N G 
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RUN VERSION 2.3 — PSR LEVEL 298 — 


GODUNOV 


000146 

000160 

000162 

000164 

000164 

000165 

000165 


IFtT.GE.TNAX. OR. CYCLE. GE.CYMAX) GO TO 9 
T = T ♦0 T 

CYCLE = CYCLE ♦ 1. 

GO TO 1 

9 CALL ENDPLTS 
CALL EXIT 
END 
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RUN VERSION 2.3 --PSR LFVFL 298 — 


000002 


000002 


000002 

000002 

000002 

000002 

000002 

000002 

000002 

000002 

000002 

000002 

000012 

000013 

000022 

000030 

000030 

000033 

000041 

000074 

000112 

000131 

000155 

000170 

000173 

000176 

000212 

000213 

000215 

000216 

000217 

000220 

000222 

000225 

000233 

000237 

000245 

000247 

000254 


SUBROUTINE INPUT 

COMMON P<50, 20,61, R(55), Z(55,205, 7START(5Q), 7ENO(50), PINTT(6> 

1, PSPOT (6) , ZOL Q ( 5 0 ) ,F POPE <50, 6) ,P A FT (5 0,6) , FORSFG1505 ,FORVEL(5 05 , 
2AFISEG(50) , AFfVEL (50) , FN120.6), PE(6), AN 120) ,UFOOF(505,UA C T(50), 
3 VNE (20) , I X Y 

4,ZAXIS (20) 

COMMON RN, RS, ZNE, ZNW, ZSE , ZSW, ESEG, WSEG, GAMMA, T, TMAX , DT, 
1 CYMAX, CPRINT, OTMIN, inT, JOT, PRCB, CYCLE, COUMP, LREAD, MWRITE 

2, LTAPE, OELR, ISPOT, VCLOLO, VOLNEW, AS, AW, AE , VNORME, VNOPMW, 

3 UNORME, UNORMW, U X ANGE, UTANGW, VFORE , VAFT, VNW, VSF,VSW, IMAX, 

4 JFORE , JAFT, I, J, PMAX, HMACHF, HMACPA, IPMAX, JPMAX 
COMMON PSPOT, ZSPOT 

DIMENSION TITLE (8) , XTITLE<8) ,YTITLE (8) 

DIMENSION ZCIRCL (55 ) , RCIRCL(55> 

INTEGER OPT 

NAMELIST /PRELIM/ IMAX, JFORE, JAFT, GAMMA, T, TMAX, CYMAX, C PP INT, COUMP 
NAMELIST/MFSH1/ RM AX, D7FORE, OZAFT, ZSPO T , DZSPO I , I SPOT , A , Z A X IS , OPT 
NAMELIST FMESH2/ R , Z , ZS T ART , ZEND , ISPOT 
NAME LIST/ST ATE 17 PIMT, PSPOT, S 
NAMELISTXSTATE2/ PIMT, PSPOT, P, PMAX,IXY 
C READ IN PROBLEM AND CYCLE NUMBER 

READ ( LPEA 0,100) CPRCP, CCTCLE 
C IS IT A RESTART OR A NEW PROBLEM* 

IF TCCYCLE .F 0 . 0.) GO TO 2 
C RESTART -READ INPUT FPOM TAPE 

1 READ (LTA«E) PROS , CYCLE 

WRITE (MWRITE, 1045 CYCLE 
104 FORMAT ( 8H CYCLE =,F7.1) 

IF (EOF, LTAPE) 14,13 

13 IF (ABS(PROB-CPROB) .GT. .01) GO TO 10 

READ (LTAPE) IMAX, JFOPE, JAFT, GAMMA, T, TMAX, CYMAX, CPRINT, 

1 ISPOT, COUMP, PMAX, DT 
READ (LTAPE) ( R ( I ) , ZS T APT( I ) , ZENO ( I ) , 1=1, IMAX) 

READ (LTAPE) ( ( Z 1 I , J) , J=1 , JFORE ) , I =1 , IMAX) 

READ (LT APE) ( < <P ( I , J , <1 , < = i , 6 ) , J= 1 , JFORE) , 1= 1 , I M AX ) 

REAR (LTAPE) (PINIT(K), PSPOT ( K > , K= 1 , 6 5 
READ (LTAPE) 

IF (CYCLE .LT. CCYCLE) GO TO 1 
14 READ (LREAO ,101) TMAX, CYMAX, CPRINT, COUMP 
GO TO 3 

C NEW PROS -READ INPUT FROM CAROS 

2 PROB = CPROB 
CYCLE = 0. 

GAMMA = 1.4 

T = 0. 

RE AO (LREAO, PRELIM) 

READ (LREAD, MESH1) 

C SET UP CLUSTERED GRID 

IF ( ( A — 1 . ) .LT. .01) R(1)=RMAX7 FLOAT (IMAX) 

IF ( ( A-l . ) .LT. .01) GO TO 65 
Rll)=(A-l.» /( ( A** IMAX) -1.) *RMAX 
65 00 66 1=2, IMAX 

R ( I) =R ( 1-1 ) *• A*(R(I-1)-R(I-2I) 

IF (I. EG. 25 R( I)=R (15 M1. + A5 
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RUN VERSION 2.3 —PSP LEVEL 298-- 


INPUT 


000261 

000264 

000266 

000270 

000271 

000300 

000312 

000326 

000346 

000354 

000361 

000364 

000366 

000370 

000370 

000401 

000402 

000403 

000404 

000405 

000407 

000411 

000413 

000414 

000420 

000421 

000426 

000427 

000431 

000432 

000434 

000452 

000461 

000465 

000470 

000470 

000471 

000471 

000501 

000507 

000513 

000514 

000527 

000532 

000536 

000536 

000546 

000547 

000554 

000564 

000574 

000606 

000617 

000620 


66 CONTINUE 

RE AO <LREAD,MESH2) 

4 DO 6 I=1,IMAX 
DO 5 J = 1 , JFORE 

IF (J .LE. JA FT) Z(T,J> = FLOAT < J) *0ZAFT 

IF (J .GT. JAFT) Z(I,J! =' r LOAT( JAETl*OZAFT «• FLO AT ( J- J AFT 1 *0 7F0P F 
IF (OPT.EO. 1) Z (I , J )=Z(I, JH-RMAX* ( ( R ( I > 7RMAX) * *2 - (R ( I 1 /RM AX) **4/2. 
1 1 

IF (OPT.EO. 2) Z (I, J)=Z (I, JM-0 .1*RMAX* ( l.-EXP(-20 . * ( PI II / RMAX ) * *2 ) » 
IF (I ,EQ. 11 ZAXIS1J) = Z(I,J> 

IF (I .EQ. 11 ZOLD (J1=Z(I,J» 

5 CONTINUE 

6 CONTINUE 

RSPOT = DZSPOT/2. 

00 17 J— 1 , J FOR E 
1/ ZOLO(J)=Z(l,J) 

PI NI T ( 3 ) = 0. 

PINIT (4) = 0. 

PSPO T ( 3 ) = 0. 

PSPOT (4) = 0. 

READ (LREAD ,S T A TE 1 1 
GM1 = GAMMA - 1. 

PINIT (51 = PINIT (1) /GM1 
PSPOT <51 = PSPOT ( 1 1 /GM 1 

PINIT ( 6) = SQRT(GAMMA*PINIT(1)/PINIT(21 ) 

IF (PSPOT ( 2 ) .LE. 0.1 GO TO 15 

PSPOT ( 61 = SORT(GAMMA*FSPOT(1>/PSPOT<2)> 

15 IF (S) 7,8,7 

7 DO 9 1=1 , IM AX 
DO 9 J = 1 , J F ORE 

IF (J .LE. JAFT 1 GO TO 11 

X = (2.*Z(I,J)-(Z(I,JFORE)*Z(I,JAFT)))/(Z(I,JFORE)-Z(I,JAFT1) 
P(I t J,l> = PINIT (1) M1.+S*X> 

IF (P(I,J,1) .LE. P MAX ) GO TO 16 
PM A X =P ( I , J , II 
IPM A X= I 
JPM AX= J 

16 CONTINUE 

P(I,J,2> = PINIM2) M1.*S*X/GAMMA1 
P ( I, J , 31 = PINIT (6) *S*X/GAMMA 
P ( I , J, 4) = 0. 

IF (I.EQ.11 GO TO 18 

ESEG=SQRT( CR(I)-R{ I — 1 ) }**2*(ZtI*J> -7(1-1, J) 1**2) 

SlNE=tR(I)-R*I-l) 1 /ESEG 
COS INE= ( Z ( I , J 1 -Z ( I - 1 , J) 1/ESEG 
GO TO 19 

18 ESEG = SORT (R ( 1 ) ** 2 ♦ ( Z ( 1 , J) - Z AX I S ( J ) ) **2> 

SINE = P(l) /ESEG 

COSINE = (Z(l, J1-ZAXIS(J))/ESEG 

19 P ( I » J , 3) = PINIK6) *S*X7GAMMA*SINE 

P ( I , J , 4 ) =-PIN I T (6) *S*X/GAMMA*COSINE 

P(I, J,5)=P(I, J,1)/GM1*P<I, J,2) *P(I,J,3)**2/2. 

P ( I , J , 6) = SQRT(GAMMA*P(I,J,11/P(I, J,21> 

GO TO 9 

11 DO 12 K=l,6 
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RUN VERSION 2.3 — PSR LEVEL 293-- 


INPUT 


000622 

12 

P ( I , J , Kl - PIN I T ( K ) 

000634 

9 

CONTINUE 

000641 

8 

READ (LREAO, STATE2I 

000644 


IF (P(l,l,2> .GT. l.E-10) GO TO 3 


C 

THE FOLLOWING DATA REFER TO THE COLLINS AND CHEN SUMP 

000650 


00 99 1=1* IMAX 

000651 


IF <R(I) .IT. 1.5663) Z(I,JFORE?= 8.437 

000657 


IF (R(I) .GT. 1.5663 .AND. R { I) .L T . 3. 1 27 ) Z ( I , JFORE > = . 85* R ( I ) ♦ 7. 2 

000675 


IF (R ( I ) . GT .3 . 127 ) 7(1, JFORE ) =10. 

000704 


Z(I, JAFT)=Z(I, JFORE 1-07 AFT 

000713 


00 20 J=l, JFORE 

000715 


IF (J .GT. JAFT .ANC. J .LT. JFORE) Z ( I , J ) =Z ( I , J A FT) * OZFO RE 
1 *FLOAT (J-JAFTI 

000736 


IF (R ( II . GT. 1.567 .AND. R(I) .LT. 3.127) GO TO 97 

000747 


P(I,J,1) = 1.126E04 

000753 


P ( I , J * 2) = .66 8 IE- 0 2 

000757 


P(I,J,3) =1667. 

000762 


P ( I * J * 41 = 0. 

000765 


GO TO 98 

000765 

97 

CONTINUE 

000765 


P(I,J,11 = 5462. 

000772 


P ( I , J, 2) = .4448E-02 

000775 


P ( I , J * 4 ) =-622. 

001000 


P(I,J,3) = 622. 

001004 

98 

Pt I, J,5)=P (I, J,l> /GMH-P (I, J,2>* (P (I, J, 31**2+ PTI, J, 4) **2) 72. 

001022 


P(I*J*6)=SQRT (GAMMA*P(I,J*1) 7P(I*J*2> ) 

001032 

20 

CONTINUE 

001035 

99 

CONTINUE 


C 

WRITE CUT INITIAL STATE 

001037 

3 

WRITE (MWRITE, 102) PROS, CYCLE 

001047 


WRITE (MWRITE, PRELIM) 

001052 


WRITE (MWRITE, MESH1 ) 

001055 


WRITE (MWRITE, STATED 

001060 


WRITE (MWRITE, MESH 2) 

001063 


RE AO (LREAO , 106 ) TITLE 

001071 


REAO(LREAO, 106) XTITLE 

001077 


READ (LREAO, 1061 YTITLF 

001105 

106 

FORMAT (SA10 ) 

001105 


Z (IMAX*2,il=0. 0 

001107 


R(TMAX+2)=0.0 

001110 


Z( IMAX+3,1) =?. 0*R( IMAX) 

001112 


R(IMAX*3)=0.0 

001113 


Z( IMAX *4, 11 =0. 0 

001114 


R(IMAX*'4)=0 .0 

001115 


Z(IMAX*5,1)=0.0 

001116 


R(IMAX*5) = R(IMAXI 

001117 


CALL XYPLOT (Z ( IMA X* 2 , 1 ) ,R ( IMAX + 2) , 4 , 1 , 0 , 0 ,0 ,1 0 . 0 , 6. 0 , XTITLE , 5 , 
1 YTITLF, 4, TITLE, 801 

001136 


RETURN 


C 

WRONG FRORLEM 

001137 

10 

WRITE (MWRITE, 1 03) PFCR,CRR09 

001147 

100 

FORMAT (2F10 .21 

001147 

101 

FORMAT (4F10.3) 

001147 

10? 

FORMAT ( 1H1 , 20 X , 15HPFCGR AM GODUNOV/ /9 X , 1 2H PROSLEM N0.,F6.1,5X, 


1 9HCYCLE N0.,F8.1> 
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INPUT 


001147 

001147 


103 FORMAT (14H WRONG PR0BLEM/16H CARO PROG NO. =,F8.3,17H TAPE PROB 
1NO . =,F8.3) 

END 
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SUBROUTINE NODES 

000002 COMMON P<50,20,6>, R(55), 7(55,20), 7START(50>, 7FND(50), PINIT16) 

1, PSP0T<6> ,ZOLD(50) ,PFOPE(50,6) ,PAFT(50,6),FORSEG<50) ,FORVFL<50) , 
2AFTSEG<50> , AFT VEL (5 0) , PN(20,6), PE<6), AN < 20 ) , U FORE ( 50 ) , UA FT t 50 ) , 

3 VNE (20) , I X Y 
4 , Z A X IS (20) 

000002 COMMON PN, RS , ZNE , ZNW, ZSE, ZSW, FSEG, WSEG, GAMMA, T, TMA X , OT, 

1 CYMAX, C°RTNT , OTMIN, IOT, JOT, PROP, CYCLE, COUMP, LREAD, MWRITE 

2, LTAPE, OELR, ISPOT, VOLOLO, VOLNEW, AS, AW, AE, VNORME, VNORMW , 
.3 UNORME, UNORMW, UTANGE, UTANGW, V FORE , VAFT, VNW, VSE,VSW, IMAX, 

4 JFORE, JAFT, I, J, PMAX , HMACHF, HMACHA, IPMAX, JPMAX 

000002 COMMON RSPOT, 7SP0T 

000002 IF (I ,GT. 1) GO TC 2 

000006 FORSEG(l) = SORT (R ( 1) **2 f t Z (1 , JFORE) -ZAXTS( JFORE) > **28 

000016 AFTSEC(l) = SQRT(R(1)**2 f ( Z < 1 , J A FT ) - 7 A X IS ( J A FT ) ) * * 2 ) 

000026 FORSEG (2 ) = SORT ( ( R (2 ) -R ( 1 ) ) **2 * (Z(2, JFORE) - Z 1 1 , J FORE ) ) * * 2 ) 

000040 A FTSE G ( 2 ) = SORT i ( R (21 -R (1 ) ) **2 f (Z(2,JAFT) - 7 ( 1 , J A FT) ) **2) 

000052 RETURN 

000052 2 IF (I .FQ. IMAX) 9E1UPN 

000055 FORSEG (1*1) =SQRT ( (R (Ifl)-R(I) ) **2 f ( Z ( If 1 , JFO RE ) -Z ( I , JFO RE > ) * *2 ) 

000072 AFTSEG(Ifl)=SORT< (R (If l)-R(I) )**2 f ( Z ( If 1 , JAF T ) -7 ( I , JAFT ) ) • * 2 ) 

000106 RETURN 

000107 ENTRY NOOES2 

00011". II = Ifl 

C COMPUTE CELL NOOE VELOCITIES 

000116 3 IF (I .EO. IMAX) GO TO 4 

000120 VNFORE=FORSEG( I)*FCRSEG (II ) 7 ( FORSE G ( I ) f FORSEG ( II ) ) * ( FORVEL (I ) f (R (I 

l)-R(I-l)) fFORVELtIl)/(P(Ifl)-R(I) )) 

0 00135 VNAFT=AFTSEG(I) *AFTSEG (111 /(AFTSEG (I) f AFISEG(Il) > *< AFTVEL (I) / (R( I) 

l-R(I-l) ) f AFTVEL (II)/ (R(Ifl)-R(I)) > 

000152 RETURN 

000153 4 VNFORE = FORSE G ( I ) * FOR VEL ( I) / ( R ( I M A X) -R l IMA X- 1 ) S 

000160 VN AFT = AFTS E G ( I ) * AF T VEL ( I) 7 < R < I M A X) - R ( IMA X- 1 ) ) 

000164 RETURN 

000165 ENTRY NOOES3 

000172 IF (J .EO. 1) GO TO 6 

000174 VNW = VNE 1 il-1 ) 

000176 VSW=VSE 

000177 VSE " VNE ( j) 

000200 GO TO 7 

000201 6 VNW = 0. 

000202 VSW = 0. 

000203 IF (I .GT. 1) VSE-VNEfl) 

000207 7 IF (J .LE. JAFT) VNE ( J ) =VN AFT *FLOA T ( J > /FLOAT ( J A F T ) 

000217 IF (J .GT. JAFT) VNE(J)=VNAFT f ( V NFORE- VN AFT) *FL OAT ( J- J A FT I /FLO AT 

1 (JFORE-JAFT) 

000233 IF ( I .EO. 1 .AND. J .LE. JAFT) VSE= A FTVEL ( 1 ) * A FTSEG ( 1 ) /R ( 1) * 

1 FLOAT(J)/FLOAT(JAFT) 

000250 IF ( I .EO. 1 .AND. J .GT. JAFT) VSE= AFT VEL ( 1 ) * AFTSE G ( 1 ) /R ( 1 ) 

1 f (FORVEL (1) *F0RSEGJ1) - A FTVEL ( 1 ) * A FTSE G ( 1))/R<1> *FL0AT(J- 

2 JAFT) 7FLOAT( JFORE- JAFT) 

C COMPUTE NORMAL VELOCITIES AT EAST AND WEST BOUNDARIES 

000275 RN=R ( I ) 

000277 ZNE-Z ( I , J) 
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NOOES 


000303 

000304 

000307 

000311 

000316 

000323 

000330 

000333 

000335 

000345 

000350 

000352 

000355 

000357 

000365 

000373 

000400 

000404 

000405 

000412 

000416 

000423 

000426 

000432 

000433 


zsw 

= ZSE 



IF 

(I 

.GT. 1) 

ZSE 

= ZOLP(J) 

ZNW 

= ZOLO ( J-11 



IF 

(CYCLE .EO. 

0. 1 

ZSE=Z ( I- 1 , J ) 

IF 

(CYCLE .EO. 

0. 1 

ZSW-ZtI-1, J-l) 

IF 

(CYCLE .EQ. 

0. 1 

7NW = Z (I ,J-1) 

IF 

( J 

.EO. 11 

ZNW = 

0 . 

IF 

( J 

.EO. 1) 

ZSW = 

0 . 

IF 

(I 

.EQ. 1 . 

AND. 

J .GT. 11 ZSW= 

IF 

(I 

• EQ. 1) 

ZSE = 

ZAXIS(J) 

RS = 

R ( I 

-11 



IF 

(I 

.EQ. 1) 

RS= 0 

• 

OELR = 

RN - RS 




ESEG =SQRT 10ELR**2«- (ZNE-ZSE) **21 
WSEG = SORT (OELR**2* (ZNW-ZSW) **2» 

VNOPME = . 5*0ELR/ESEG* ( VNE ( J) *VSE) 

VNORMW = . 5*OELR7WSEG* (VNW+VSW) 

RETURN 

ENTRY NOOES4 

COMPUTE NEW NODE POSITIONS, I.E. MOVE MESH 
ZOL D ( J 1 - Z < I, J ) 

Z(I,J) = ZII,J) ♦ OT*VNE(J) 

IF T I .GT. 1) RETURN 

ZA XI S < J 1 = zaxisu) ♦ DT*VSE 

RETURN 

END 
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000002 

000002 

000003 

000011 

000013 

000014 

000016 

000020 

000027 

000054 

000066 


000114 

000125 

000133 

000135 

000137 

000142 

000142 

000143 

000153 

000163 

000166 

000170 

000170 

000204 

000204 

000205 

000206 

000214 

000215 

000217 

000217 

000221 

000222 


SUBROUTINE SHOCK 

C SHOCK SOLVES THE FULL RIMANN PROBLEM AT THE FORE AND 

C AFT SHOCKS TO DETERMINE THE SHOCK VELOCITY AND FLUID 

C PROPERTIES BEHIND THE SHOCK 

COMMON /SH0CK/P1,RH01,U1, VI , P4 , RH04 , U4 , 'V4 , G , V SHOCK , P3 , RH03 , U3 , M 

GAMMA=G 

C IS THERE REALLY A SHOCK WAVE 

IF ( ABS(1.-P1/P4) .LT. l.E-10) GO TO V 
C START EY GUESSING PRESSURE ACROSS CONTACT SURFACE 

PCS =(Pl*P4*/2. 

C COMPUTE FLUX DENSITY OF W4 (SHOCK WAVE! 

ITER -0 

GMi = GAMMA-1. 

GP1 = GAMMA H. 

1 EM4 =SQRT(RH04/2.* < C-Pl *PCS *-GMl *P4 ) ) 

C COMPUTE FLUX DENSITY OF W1 (RAREFACTION WAVE* 

EMI = GMl/2./GAMMA*SQRT (G AMM A*P1* R HOI > * ( 1 . -PCS/PI ) / ( 1 . - ( PCS/ PI ) * * 
KGM1/2. /GAMMA) ) 

IF (P1.LT.P4I EM4=SCRT TPH01/2.* (GP1*PCS+GM1*P1 > > 

IF (P1.LT.P4) EMl = GMl/2./GAMMA*SQRT<GAMMA*P4*RH04)Ml.-PCS/P4) / 

1 (1.- (PCS/P4) ** (GM1/2./GAMMAM 
C CHECK GUESS FOR PCS 

PCSP = (EM1*P4 + EM4*P1 + EM1*EM4* (U1-U4* > t CEM14-EM4 1 

IF (ABS(1.- PCS /PCSP* . LE. l.E-03) GO TO 2 

ITER = ITER 1 

IF ( ITEP . GE. 50) GC TC 3 

PCS = (PCSP*PCS» /2. 

C HAVE NOT CONVERGED, GUESS AGAIN AND REPEAT CYCLE 

GO TO 1 

C CONVERGENCE, COMPUTE JUMP CONDITIONS ANO SHOCK VELOCT 

2 P3 = PCSP 

RH03 = RH04* <GP1*P3«-GM1*P4) / ( GP1 * P 4+ GM 1*P3 5 
U3 =(P1 -P4«-EM4*U4 + EM1*U1) /(EM1»EM4) 

VSHOCK = U4 EM4/RH04 

V3 = V4 

RETURN 

3 WRITE (M , 1 0 0) PI, P4, PCS, PCSP 

100 FORMAT (44H SHOCK HAS NOT CONVERGED AFTER 50 ITERATIONS/ 17H PI, 04 
1, PCS, PCSP = ,4E 15.4) 

CALL ENDPLTS 
CALL EXIT 

4 VSHOCK = U44-SORT (G*F4/RH04) 

RH03=RH04 
P3 = P4 
U3=U4 
V3 = V4 
RETURN 
END 
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000002 


000002 


000002 

000002 

000002 

000004 

000006 

000011 

000013 

000015 

000016 

000017 

000023 

000026 

000030 

000032 

000036 

000040 

000044 

000044 

000046 

000053 

000065 

000077 

000111 

000121 

000125 

000130 

000133 

000142 

000175 
00021 3 
000232 
000256 
000271 
000303 
000311 
000312 
000317 
000325 
000332 
000334 
000335 
000342 


SUBROUTINE OUTPUT 

COMMON P(50, 20,6) , fi{55>, Z<55,20>, ZSTART(50), ZEND(50>, PINIT16) 

1, PSPO I (6) , ZOLOC50 ) ,PFCRE 150,6) ,Pft FT ( 50 , 6» , FORSE G (5 0) ,FORVEL(5 0) , 
2AFTSEG150) , AFTVEL(50) , PN(20,6)» PE<6), A N ( 20 ) , U FORE( 50) , UAET < 50 ) , 
3 VNE (20) ,IXY 

4,7AXIS(20) 

COMMON RN, RS, ZNE , ZNW, ZSE, ZSW, ESEG, WSEG, GAMMA, T, TMAX , DT, 
1 CYMAX, CPRINT , OTHIN, IOT, JOT, PROB, CYCLE, CDUMP, LREAO, MWRITE 

2, LTAPE, OELR, ISPOT, VOLOLD, VOLNEW, AS, AW, AE , VNORME , VNORMW, 

3 UNORME, UNORMW, UTANGE, UTANGW, VFORE, VAFT, VNW, VSE,VSW, IMAX, 

4 JFORE , JAFT, I, J, PMAX, HMACHF , HMACHA , IPMAX, JPMAX 

COMMON RSPOT, ZSPOT 

DIMENSION ZPRIME (55 ,2 ) ,PPRIME (55) 

MIMAX=-IMAX-1 
CPLOT = 2. *CPRINT 
IF (AMOD(CYCLE, CPLOT) ) 25,16,25 
16 ZPRIME (1,11 = ZAXIS (JFORF) 

ZPRIME (1,2) = ZAXIS (JAFT) 

RPRIME (1 ) = 0 .0 
00 14 1=1, IMAX 

Z°R I ME (1 + 1,1 ) =Z (I , JFORE) 

ZPRIME (IH, 2 ) =Z ( I, JAFT) 

RPRIME (IH»=R(I) 

14 CONTINUE 

CALL XYPLOT (ZPRIME (1,1 ), RPRIME 1 1 ), M IMAX , 1 , 0 , 1 ) 

IF (JAFT.EQ.l) GO TC 33 

CALL XYPLOT (ZPRIME ( 1,2 > ,RPRI ME < 1 ) , M I MA X , 1 , 0 , 3 ) 

33 CONTINUE 

IPLOT = IPLOT +1 
WRITE (MWRITE, 115) CYCLE 
25 WRITE (MWRITE, 100! FROR, CYCLE, T 
WRITE (MWRITE, 101) TT, IOT, JOT 
WRITE (MWRITE, 102) F MAX, IPMAX, J PM AX 
WRITE (MWRITE, 103) Z A X IS < J FORE ) , ZAXIS(JAFT) 

IF (AMOD(CYCLE ,COUMF) 1 1,3,1 

1 IF (T-TMAX) 2,3,3 

2 IF (CYCLE-CYMAX) 4,3,3 

3 WRITE (LTAPE) PROB , CYCLE 

WRITE (LTAPE) IMAX, JFORF, JAFT, GAMMA, T, TMAX, CYMAX, CPRINT, 

, 1 ISPOT, COUM 3 , PMAX, OT 

WRITE (L TA PEI (P ( I ) ,ZS T ART ( I ) , ZENO < I ) , 1 = 1, IMAX) 

WRITE (LTAPE) ( ( Z t I , J ) , J=1 , JFORE) , 1 = 1 , IMA X ) 

WRITE (LTAPE) (((P(I,J,K>, K=1 , 6) , J=1 , JFORE ) , I =1 , IMAX I 
WRITE (LTAPE) (PINTT(i<>, PS°OT (<), <=1,6> 

WRITE (LTAPE) ( Z OL C ( J ) , J- 1 , J F ORE) 

WRITE (MWRITE, 111) CYCLF 

4 RETURN 
ENTRY OUTLNG 

IF (IXY.LE.O) WRITE (MWRITE, 104) 

IF (IXY.GT.OI WRITE (MWPITE, 114) 

DO 11 1=1, IMAX 

00 10 J=i, JFORE 

IF (J .EO. 1) WRITE (MWRITE, 110) 

WRITE (MWRITE. , 105) I, J, R(I), 7(1, J), P(I,J,1), P(I,J,2), P(I,J,3 
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OUTPIJT 


1* , P< I,J,4) , P < I, J , 5) 

000401 10 CONTINUE 

000404 11 CONTINUE 

000406 WRITE (MWRITF,117) 

000412 WRITE (MWRITF,116> C J, Z AXIS! J) , J=1 , JFORF) 

000426 M=MWRI TF 

000430 WRITE CM, 1061 

000433 WRITE CM, 1071 

000437 DO 12 1 = 1 , IMAX 

000441 12 WRITE CM, 109) I , FOR VEL ( I 1 , HM ACHF , { PFOPE J T , K 1 , K=1 , 61 

000465 WRITE (M,103) 

000470 WRITE CM, 107) 

000474 00 13 1=1, IMAX 

000476 13 WRITE (M,109) I , A F T VEL ( I 1 , HM ACHfl , < RAFT < I , K > , K = 1 , 5 ) 

000522 100 FORMAT { 10 X , 3H FRO 0LFM , F7 . 2 , 6H CYCLE, r B. 1 , 4H T =,F12.4t 

000522 101 FORMAT (/5H OT =,E12.4,12H AT CELL I =,I3,4H J =,I3> 

000522 102 FORMAT C/7H PMAX =,F12.5,12H AT CELL I =,I3,4H J =,131 

000522 103 FORMAT <10X,21H FORE SHOCK IS AT Z =,El2.4,20H AFT SHOCK IS AT 7 = 

1 , E12 . 4 1 

000522 104 FORMAT (4 X , 1HI , 4X , 1 H J , 9 X , 1HR, 14X,1HZ,14X,1HP,12X, 3 HRHO , 14 X , 1 HU , 14 X, 

11H V , 9 X , 6HENERG Y 1 

000522 105 FORMAT (215, 7(3X,E12.4)1 

000522 106 FORMAT <1H1,20X,22H FORE SHOCK PROPERTIES) 

000522 107 FORMAT ( 7X , 1H 1 , 3 X , 6 H VS HOCK , 7 X , 6HMS HOC K , 12X , 1H® , 1 2 X, 3HRHO , 14X , 1 HU , 

114X,1HV,10X,6HENERGY) 

000522 108 FORMAT (20X.21H AFT SHOCK PROPERTIES) 

000522 109 FORMAT ( 5X , 15 , 7 T3X , E12 . 4) ) 

000522 110 FORMAT (IX,/) 

000522 111 FORMAT C//19H TAPE TUMP ON CYCLE, F7.17) 

000522 112 FORMAT C/7H PMIN =,E12.4,12H AT CELL I =,I3,4H J =,I3> 

0 0 0522 114 FORMAT (4 X , 1 HI , 4 X , 1 H J , 9 X , 1HY, 14X , 1 H X , 14 X , 1HP ,1 2 X , 3HRHO , 14X , 1 HU , 14 X, 

1 1 H V, 9X , 6HE NERG Y) 

000522 115 FORMAT (1H1,5X,37H FCRE AND AFT SHOCKS PLOTTED AT CYCLE, F3.1) 

000522 116 FORMAT ( 4 ( 15, 1 X , E 1 2 .4 ) ) 

000522 117 FORMAT (12H ZAXIS ARRAY) 

000522 RETURN 

000522 ENO 
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000012 

000020 

000030 

000030 


SUBROUTINE EOS < G , P , f HO , ZOO T , ROOT , E NERGY , SOUND > 
P=(G-1.)*(ENERGY - RH0M700T**2 ♦ ROOT * *2 ) 72. ) 
SOUNO = SORT(G*PZRHC> 

RETURN 

END 
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C 

C 

000002 


000002 


000002 

000002 

C 

000002 
000004 
000006 2 
000015 
000015 5 

000017 6 

C 

000022 10 
000025 
000030 
000034 

C 

000037 

000041 

000042 12 

000045 

000046 15 

000051 

000054 


SUBROUTINE FLUXES 

COMPUTES MASS, RADIAL MOMENTUM, AXIAL MOMENTUM, AND ENERGY 
FLUXES ON THE NORTH, SOUTH, EAST, AND WEST BOUNDARIES 
COMMON P <50 *20 , 6) , R<55l, Z<55,20», ZSTART (50 ) , ZEND<50», PINITC6) 

1, PSPOT (6) ,ZOLD<50) ,PFORE<50,6» , PA FT < 50 , 6 > , FORSE G (50J ,FORVEL(50) , 
2AFTSEGT50) ,AFTVEL<50) , PN<20,6), PE<6>, AN < 20 ) , U FORE < 50 ) , UAFT < 50 > , 
3 VNE ( 2 0 1 , I XY 

4,ZAXIS<20) 

COMMON RN, RS, ZNE , ZNW, ZSE, ZSW, ESEG, WSEG, GAMMA, T, TMAX , DT, 
1 CYMAX, OPR IN T , DTMIN, IDT, JOT, PROB, CYCLE, CDUMP, LREAD, MWRITE 

2, LTAPE, DELR, ISPOT, VOLOLD, VOLNEW, AS, AW, AE , VNORME, VNOPMW , 

3 UNORME, UNORMW , UTANGE, UTANGW, VFORE , VAFT, VNW, VSE,VSW, IMAX, 

4 JFORE , JAFT, I, J, PMAX, HMACHF, HMACHA , IPMAX, JPMAX 
COMMON RSPOT, ZSPOT 

COMMON/FLUX/ FLUXN<20,5>, FLUXS (5) , FLUXW <5 ) , FLUXE < 5 > 

SOUTH FCUNPAR Y 
IF (I . EQ . 1) GO TO 5 
00 2 K=2 , 5 

FLUXS(<)=-FlUXN(J,<) 

GO TO 10 
DO 6 K = 2 , 5 
FLUXS(K)=0. 

NORTH BOUNDARY 
FLUXN< J,2) = -PN ( J, 2) *PNtJ,4) 

FLUXN < J, 3>=-PN ( J, 21 *PN ( J,4)*PN ( J, 3 ) 

FLUXN(J,4)=-PN(j,21*PN(J, 41**2 - PN(J,1> 

FLUXN <J,5) = -(PN(J, 5 M-PN (J, 1> ) *PNt J,4) 

WEST BOUNDARY 
IF (J . EO . 1) GO TO 15 
00 12 <=2,5 
FLUXW(K)=-FLUXE<K> 

GO TO 16 

FLUXW (2) =P< I, 1 ,2) *F (I, 1 ,3) 

FLUXW <3>=P < 1, 1 , 2) *P <1 , 1,3) ** 2 + P < 1 , 1 , 1 ) 

FLUXW(4)=P(I,1,2) * D (1,1,4) *P(I,1,3) 


000056 
000061 16 
000064 
000072 
000103 

oooio r 
000110 


FLUXW(5)=(P< 
FLUXE (2) = - 
FLUXE (3) = - 
FLUXE (4 ) = - 
FLUXE <5) = - 
RETURN 
END 


1,1,5) ♦ P1I,1,1))*P(T,1,3) 

PE <21 * < LNORME- VNORME) 

PE (2) *PE <31* (UNORME-VNORMF > - PE < 1 ) * DELR/ ESEG 
°E (2) *PE (4) *<UNO p ME-VNORME> + PE< 1 >*< ZNE- ZSE ) /ESEG 
P E (5 ) * < LNORME- VNORME 1 - PF T 1 ) ’UNORME 
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SUBROUTINE DTCALC 

C COMPUTE TIME STEP FROM COUPANT CONDITION 

0000 02 COMMON P(50,20,6), F (55) , Z 1 55 , 201 , ZSTART(50>, ZEND(50), PINITT6) 

1, PSPOT ( 6) ,ZOL 0(501 ,PFORE(50,6) ,P A FT (5 0,6) ,FORSEG(50) ,FORVEL(50> , 
2AFTSEG (50) , AFTVEH5G) , PN(20,6), PE(6>, AN ( 20 ) , UFORE ( 50 ) , UAET ( 50 ) , 
3 VNE (20) , IXY 

4,ZAXIS(20) 

0000 02 COMMON PN, PS, ZNE , ZNW, ZSE, ZSW, ESEG, WSEG, GAMMA, T, TMA X , DT, 

1 CYMAX, CPRINT, DTMIN, IDT, JOT, PRQB, CYCLE, CDUMP, LREAD, MWRITE 

2, LTAPE, DELR , ISPOT, VOLOLO, VOLNEW, AS, AW, AE , VNORME , VNORMW, 

3 UNORME, UNORMW, UTANGE, UTANGW, VFORE , VAFT, VNW, VSE,VSW, IMAX, 

4 JFORE, JAFT, I, J, PMAX, HMACHF, HMACHA , I»MAX, JPMAX 

000002 COMMON RSPOT , ZSPOT 

000002 NAMELIST/COT/ 0 T1 , I , J , DTR, DT Z 

C RADIAL VELOCITY TIME STEP 

000002 OR=D£LR 

000004 RSIG =AMAX1 ((P ( I, J, 6) + P (I, J,4> ) , (P (I, J,6) -P1I , J, 4) ) ) 

00001Z DTR = DR/RSIG 

C AXIAL VELOCITY TIME STE D 

000021 DZ = AMIN1( (ZNE-ZNW), (ZSE-ZSW)) 

000030 ZSIG ■= AMAX1( (P( I, J .6H-ARS (UNORME) «-ABS (VNORME) ) , (P( I, J,6> +ADS t 

1 UNORMW) ♦ABS(VNORMW)) ) 

000050 DTZ - OZ'ZSIG 

C TIME STEP FOR CELL I, J 

000052 DTI = OTR’DTZ/ (DTRf ETZ) 

000055 IF (DTI .LE. 0.) GC TO 6 

000056 IF (I.EO.l .AND. J.EQ.l) OTMIN=DT 1 

000066 IF (DTI , GT . DTMIN) GO TO 5 

000072 DTMIN * OT1 

000072 IOT = I 

0000 74 JO T = J 

000075 5 RETURN 

000076 6 WRITE (MWRITE, COT) 

000101 CALL ENDPLTS 

000102 CALL EXIT 

000103 END 
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RUN VERSION 2.3 --PSR LEVEL 29°.-- 


C 

c 

000002 


000002 


000002 

000002 

000002 


000002 

000004 

000006 

000011 

000015 

000022 

000026 

000033 

000037 

000037 

000041 

000044 

000047 

000051 

000054 

000055 

C 

000063 

000102 

000112 

000114 

000114 


SUBROUTINE NEWELO 

COMPUTES CEIL SURFACE AREAS AND VOLUMES 

AND UPDATES FLOW VARIABLES WITH THE DE FFER.ENCE EONS. 

COMMON P(50,20,6), R<55), 7(55,20), ZSTART(50), TEND (50) , P T N T T ( 6 ) 

1, PSPOT (6) , ZOLDK5 0 ) ,PE0RE(5D ,6) ,PAET (50,6) ,EORSEG(50) ,FORVEL (50) , 
2AETSEG(50) ,A C TVEL (50) , PN(20,6>, PE(6), AN ( 20 ) , U FORE ( 50) , UA FT ( 50 ) , 
3 VNE (20) ,IXY 

4.ZAXIS (20) 

COMMON RN, RS, ZNE , ZNW, ZSE, ZSW, ESEG , WSEG, GAMMA, T, TMA X , Of, 
1 CYMAX, CP R TNT , DTMIN, IDT , JDT , PROB, CYCLE, CDUMP, LRE A 0 , MWRITE 

2, LTAPE, DELR, ISPOT, VOLOLD, VOLNEW, AS, AW, AE , VNORME, VNORMW , 

3 UNOPME , U NORM W , UTANGE, UTANGW, V FORE , VAFT, VNW, VSE,VSW, IMAX , 

4 JFORE, JAFT, I, J, PMflx, HMACHF, HMACHA , IPMAX, JPMAX 
COMMON PSPOT, ZSPOT 

COMMON/FLUX/ FLUXN(20,5), FLUXS(5), FLUXW ( 5 ) , F LU XE X 5) 

NAMEL 1ST 7FLUX / I, J , fi A TE , OM ASSO , VOL NE W , VOLOLO, FLUXN,FLUXS , FLUXW, 

1 FLUXE, PN, PE, UNORME, VNORME, ZNE , Z SF , ZN W, ZS W, A PL A NE , AN , A S , RN , PS 
2, AW, AE, DELR, ESEG 

SURFACE A pf AS 

PI = 3.14159 
IF (IXY.GT.O) GO TO 1 
IF (I ,EO. 1) A S= 0 • 

IF (I .GT. 1) A S = A N f J ) 

AN ( J) =PI*2. *RN* (ZNE -ZNW) 

IF (J .GT. II AW = AE 
IF (J .EO. 1) AW = PI*<RN+RS)*OELR 
AE = PI* (R N *-RS ) *E S E G 
GO TO 2 

1 AS= AN ( j I 

ANt J)=ZNE-7NW 
IF (I.EO.l) A S = A N ( J ) 

AW = AE 

IF (J.EO.l) A W = DELE 
AE=ESEG 

2 aplane= o • 5 *oelr* ( Z NE-ZNW* ZSE-ZSW ) 

OLD CEIL VOLUME 
IF (IXY.LE.O) \JOL= 

1 PI/3 . *DELR* ( ( 2 . *KS»PN) * (ZSE-ZSH) * (2. *RNf RSI * (ZNE-ZNW ) > 

IF (IXY.GT.O) VOL- (EN-PS) MZNE-ZNWfZSE-ZSW) 72. 

VOLOLD= VOL 
RETURN 

ENTRY NUELOW 

NEW CELL VOLUME 


000121 


ZNNE = Z T I , J) 


000125 

ZNSW-Z (1-1, J- 

11 

000130 

ZNSE=Z ( 1-1 , J) 


000133 

ZNNW=Z (I, J-i) 


000136 

IF 

(I .EQ. 1) 

ZNSW = ZAXIS(J-l) 

000141 

IF 

(I .EO. 11 

ZNSE = ZAXISTJJ 

000144 

IF 

(J .EQ. 1) 

ZNNW =0 « 

000147 

IF 

(J .EQ. 1) 

ZNSW = C . 

000151 

IF 

( IXY.LE.O) 

VOL = 


1 

PI/3.*D£LR* ( (2. *RS*RN) MZNSE 

000171 

IF 

(IXY.GT.O) 

VOL= ( EN-RS) * IZNNE- 
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RUN VERSION 2. 3 * - PER LEVEL 293-- 


NE WEL 0 


000201 

C 

000203 
000204 
000211 
000212 
000229 
000227 
000231 
000240 10 

000261 

000306 20 

000310 

000317 

000337 

000367 

000360 


VOLNEW = VOL 

SOLVE CONSERVATION EONS FOR NEW FLUIO PROPERTIES 
CALL FLUXES 

OMASSO = P(I, J,2> ’VCLOlO 
00 20 K=2 , 5 

RATE = AN( J)*FLIIXN(J,iO+AS*RLUXSt<)+AE*FLUXE(K» + AW*FLUXW(K) 

IF (X .NE. 41 GO TO 10 

IF ( I X Y .GT. 01 GO TO 10 

p.AIE = RATE «■ 2.»PI*A°LANE*P(I, J, 1 > 

IF (K.F0.2 .OR. < . E C. 5 ) P ( I , J , < ) = t VCLOL D*P ( I , J , < ) *R AT E *DT ) / VOLNE W 
IF (K.E0.3 .OR. K.EC.4) P ( I , J , K! = l QMASS0*° 1 1 , J , K ) «-R A TE*OT ) f° < I , j , 
1 2 ) /VOLNEW 
CONTINUE 

IF <P(I,J,5) .LT. 0.1 WPITE <6, FLUX) 

IF (I .EO. 1 .AND. J .EO. 7 .AND. CYCLE .LE. 2.) WRITE (6, FLUX) 

IF (I .EO. 1 .AND. J .EO. 1 .ANP. CYCLE .LT. 2.) WRITE 16, FLUX) 

RETURN 

END 



RUN VERSION 2.3 — PSR LEVEL 298-- 



C 

SUBROUTINE RIO 

RIO EVALUATES THE PLOW VARIABLES ON THE EAST 


c 

AND NORTH CELL BOUNDARIES BY SOLVING A RIEMANN PROBLE 

000002 


COMMON P (50 ,20 ,6) , F<55), Z<55,20), ZST ART (50 ) , ZEND(50), PINIT(6> 

000002 


1, PSPOT (6) ,70LD (50 ) , F FORE ( 50 , 6 > , PA FT t 5 0 , 6 ) , FOR SE G 150 ) ,FORVEL (5 0) , 
2AFTSEG (50) » AFTVEH5G) , PN(20,6)> PE(6), AN ( 20 ) , UFORE ( 50 ) , UAFT ( 50 ) , 
3 VNE (20) , I X Y 
4,7AXIS (20) 

COMMON RN, RS, ZNE, ZNW, ZSE, ZSW, ESEG, WSEG, GAMMA, T, TM AX , OT, 

000002 


1 CYMAX, CPRINT , DTM IN , IOT , JOT, PROB, CYCLE, COUMP, LREAD, MWRITE 
2, LT APE, OELR, ISPOT, VOLOLD, VOLNEM, AS, AH, AE , VNORME , VNOPMW , 

3 UNORME, UNORMW, UTANGE, UTANGW, VFORF , VAFT, VNW, VSE,VSW, IMAX, 

4 JFORE, JAFT, I, J, PMAX, HMACHF, HMACHA, IPMAX, JPMAX 

COMMON RSPOT, 2SP0T 

000002 


C0MM0N/RIEMANN/P1,U1,PK01, V1,C1,P4 ,U4 ,RH04, V4,U,C4, G,PB,UB,RHOB , 

000002 


1V9 

G=GAMMA 


c 

NORTH BOUNDARY 

000004 


gmi-gamma-i. 

000006 


IF (I .EQ. IMAX) GO TO 10 

000010 


PI = P(I,J,1) 

000013 


U1 = P ( I , J , 4) 

000016 


RH01=P(I,J,2) 

000021 


V1=P(I, J,3» 

0 000 2** 


Cl = P(I,J,6) 

000022 


U = 0. 

000030 


P4-P(I*1 t J i U 

000033 


U4 = P ( I *• 1 , J , 4) 

000036 


RH04=P (Ifl, J, 2) 

000041 


V4 = P ( I*1 » J » 3) 

000044 


C4=P(I f 1 , J , 6) 

000042 


CALL RIEMANN 

000050 


PN(J,1)=PB 

000052 


PN ( J , 2 ) = RHOB 

000054 


PN(J,3)~VB 

000055 


PN( J,4)=UB 

000052 


PN ( J, 5) =P3/GMl +■ RHCS/2.»(U3»*? + V3 ,I ‘*2) 

000066 


GO TO 15 

000062 

10 

00 11 <=1,5 

000071 

11 

PN ( J , K ) = P(IMAX,J,M 


c 

EAST "CUNDARY 

000105 

15 

IF (J .EO. JAFT) GC TO 25 

000107 


IF (J .EQ. J e ORE ) GC TO 30 

000111 


SINE = OELR/ESEG 

000113 


COSINE = (ZNE-ZSE) /ESEG 

000115 


PI = P(I,J,1) 

000121 


RH01 = P (I, J , 2) 

000124 


P4=P(I,J+t,l) 

000127 


RH04 = P (I , J*l,2) 

000132 


U1 = P(I, J,3) *SINE - P(I,J,4)*COSINE 

000141 


U4 = P(I, J+1,3)*SINE - P(I,JH,4)*COSINE 

000147 


VI = F II , J,3> *COSINE * P( I,J,4)*SINE 

000155 


V4 = P (I , J ♦ 1 , 3 1 * COS INF * PfI , J+1,4) *SINE 

000162 


Cl = F(I,J,6) 
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RUN VERSION 2.3 — PSR LEVEL 293 


RIO 


000165 


C4 = P(I,JH,6) 

000166 


U = VNORME 

0001/0 


CALL riemann 

000171 


UNORMW = UNORME 

000173 


VNORMW = VNORME 

000174 


IF (J .EO. 1) UNORMW = P(I,J 

000202 


IF (J .Et2. 1) VNORMW = 0. 

000204 


UNORME = Ufl 

000206 


PE T 1 ) = PB 

00020/ 


PE ( 2 ) = RHOB 

000211 


PE<3» - UB*SINE *■ VE*COSINE 

000214 


PE <4 ) = VB * SINE - UB*COSINE 

000217 


PE 1 5) = PB/GM1 ^ RHCB/2.*(UB 

000225 


RETURN 

000226 

25 

00 26 K= 1 i 5 

000230 

26 

PE ( K > = PAFT(I,K» 

000240 


UNORMW=UNORME 

000241 


VNORMW” VNORME 

000243 


UNORME = U AFT t I i 

000244 


RETURN 

000245 

30 

00 31 <=1,5 

000247 


PE ( K ) = PFORE ( I , K) 

000255 

31 

CONTINUE 

000257 


VNORMW=VNORME 

000260 


UNORMW=UNORME 

000262 


UNORME = UFORE (I) 

000263 

5 

RE TURN 

000264 


END 


VB**2T 
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PUN VERSION 2.3 - -PSP LEVEL 293- 


0000 02 


000002 


000002 
000002 
000002 
000005 
0000 0 2 
000011 

000013 

000017 

000022 

000025 

000030 

000032 

000034 

000036 

000045 

000045 

000047 

000054 

000064 

000072 

000074 

000076 

000100 

000104 

000105 

000107 

000112 

000114 

000116 

000117 

000123 

000126 

000137 

000142 

000146 

000146 

000147 

000152 

000162 


SUBROUTINE SHIO 

COMMON P (50 ,20 ,6) , M55), 2(55,20), 7START<50), ZEND(50>, PINIT<6> 

1, PS POT (6) , 70LO <5 0 ) ,PFO'PE(50,6> ,P A FT < 5 0 , 6 ) , FORSE G (5 0 > ,FORVFL<50> , 
2AFTSEG dO) , AFTVEL150) , PN(20,6>, RE(6>, AN T 20 ) , U FORE ( 5 0 ) , UA F T 1 50 5 , 
3 VNE (20) , IXY 

4,7AXIS(20) 

COMMON RN, RS, ZNE , ZNK , ZSE, ZSW, FSEG , WSEG, GAMMA, T, TMA X , OT, 
1 CYMAX, CPRINT, DTM1N, IOT , JDT , PROP , CYCLE, COUMP, LREAD, MWRITE 

2, LTAPE, OELR, ISPCT, VOLOLO, VOLNEW, AS, AW, AE , VNOPME , VNOPMW, 

3 UNORME, UNORMW , U T ANGE , UTA NGW , VFOP.E, VAFT, VNW, VSE,VSW, IMAX, 

4 JFORF, JAFT, I, J, PMAX, HMACHF, HMACHA, I °M A X , JPMAX 
COMMON RSPOT, ZSPOT 

COMMON XSHOCK^Pl ,RHO 1,U1 , VI ,P4 , RH04 , U4 , V4, G , VSHOCK , P3 , RH03 , U3 , » 
IF (I ,EO. IMAX) RETURN 
G=GAMP A 
II=I«-i 

IF (I .FO. 1) 11=1 

SET UP AFT SHOC< INPUT FOR SHOCK 

1 P4 = P(II, JAFT*1,11 
PI = P(II, JAFT ,1) 

RH04 = P(I1,JAFT+1,2) 

RHOl = P(II,JAFT,2) 

OELR = P(II)-R(II-1) 

STHETA = OELR/AFTSEO(II) 

IF (II .EQ. 1) GO TO 2 

CTHETA = (Z(II,JAFT)- 7 (II-l , JAFT ) > 7AF TSEG ( II ) 

GO TO 3 

2 STHETA = R(l) f AFTSEG(l) 

CTHETA = (Zll, JAFT) -ZAXISf JAFT) )/AFTSFG(l> 

3 U1 = P(II,JAFT ,3) *STHETA- P( I I , JAF T , 4 ) "CTHETA 

U4 = P (II, JAFT*1,3) "STHETA- P ( 1 1 , J AFT+ 1 , 4 ) *CT HET A 
U AFT ( I I ) = U1 
Cl = P (I I, JAFT , 6) 

C4 = P(II, JAFT»1,6) 

V4 = P (I I, JAFT + 1 , 3 ) *C THET A ♦ P (I I , JAFT*1 ,4)*S THETA 
CALL SHOCK 
V3 = V4 

PROCESS AFTSHOCK OUTPUT FROM SHOCK 
IF (PROD .E0.19.) VSHOCK=0. 

AFTVEL(II)=VSHOCK 
PAFT (11,1) =P3 
PAFT ( 1 1 , 2) =PHO 3 

PAFT (11,3) =U3*STHETA * V3*CTHETA 
PAFT (11,4) =V3*STHE TA-U3*C THETA 

PAFT (11,5) = P3/(GAMMA-1.) ♦ RH037 2. * ( U3*»2*V3 ** 2) 
HMACHA=(VSH0CK-U4) / C4 

IF J=l, DO THE SECONO TUBE ALSO 
IF (II ,GT. 1) GO TC 4 
II = 2 
GO TO 1 

4 IF (I .EQ. 1) II = 1 

FORE SHOCK ANGLE 

7 CTHETF = (ZUI,JFORE)-Z(II-l, JFORF ) ) KFORSEG (I I ) 

STHETF = 0ELR7 FORSE G ( 1 1 ) 
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RUN VERSION 2.3 - - p SR LEVEL 298-- 


S H 10 


000164 

000170 

000176 

000202 

000205 

000214 

000216 

000223 

000226 

000230 

000231 

000232 

000234 

000256 

000264 

000271 

000277 

000300 

000302 

000304 

000305 

000307 

000311 

000317 

000321 

000324 

000325 

000326 


IF (II .EQ. 1) ST HE I F = 9 ( 1 1 t FO RSE G ( 1 ) 

IF (II .EQ. 1) CTHETF = ( 7 Cl , JFORE 1 - Z A X IS ( JFO R El ) 'FOPSEG ( 1 ) 
SET UP FORE SHOCK INPUT TO SHOCK 
PI = P ( I I , JFORE, U 
RHOl = P(II, JFORE, 21 

U1 = P(II, JFORE, 3) 'STHETF - P ( 1 1 , JFOP.E , 4) »CTHE TF 
UFORE ( IU = U1 

VI = P<II, JFORE, 4) *STHETF ♦ P { 1 1 , J FORE , 3) *CTH E TF 
Cl = P(II, JFORE, 61 

ARE KE IN COLO SPOT OR UNDISTURBED REGICNa 

5 P4-PINIT (1 ) 

6 U4 = 0. 

V 4 = 0. 

IF (RSPOT.LT. l.E-41 GO TO B 

RH04 = PINIT(2) * (1 . »PSP0T(2> *EXP(-3. * t (Z (I, JFORE 1-ZSPOT1 **2 
1 + R ( 1 1 ‘*21/RSF0T**2> - PSPOT ( 2 1 *E XP < -3. >) 

8 RADSOR = (Z(I, JFORE 1-7SPOT1 **2 «-R(I)**2 
IF (PAOSQR .GT. p SPCT**21 RH04=PINTT ( 21 
C4=PIN IT (61 *SQ9T (PINIT (21/RH041 

r A | i c unpi/ 

PROCESS FORE SHOCK OUTPUT 
FORVEL (II) =VSHOCK 
PFORE ( I I , 1 1 = P3 
PF ORE (11,21 = RH03 
PFORE (11,3) = U3*STEETF 
PFORE (IT ,4) =-U3*CTHETF 

PFORE ( II ,5) = P37 (GAHMA-1.) ♦ PH03K 2. *U3* *2 

HMACHF = VSH0CKKC4 

IF (II .GT. 1) RETURN 

II = 2 

GO TO 7 

END 
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RUN VERSION 2.3 — PSR LEVEL 296 


C 

000002 


000002 

000003 

000007 

000014 

000025 

C 

000036 
000046 
000057 
000062 
000064 
000066 
000067 
000070 
000072 
000074 
000074 5 

000075 
000077 
000100 
000102 
000103 10 

000104 
000106 
000114 
000116 
000117 15 

000120 
000122 
000130 
000132 
000133 


SUBROUTINE RIEMANN 

SOLVE WEArf! WAVE RIEMANN PROBLEM 

C0MM0N/RIEMANN/P1,U1,RH01, Vi , Cl , P4 , U4 , RH04, V4 , U , C4, G ,PB , UB, RHOB , 
1VB 

COMPUTE PRESSURE AND VELOCITY ACROSS CONTACT SURFACE 

GAMMA=G 

EMI = SORT (GAMMA*P1 *RH01> 

EM4 = SORT <GAMMA*P4«RH041 

UCS =(Pi - P4 ♦ EM1*L1«-EM4*U4)/(EM1«-EM4> 

PCS = <EM1*EM4MU1-1'4) » P4*EM1 ♦ Pi*EM4) / 1 EM 1 +E M4T 
COMPUTE WAVE SPEEDS 


VW1 

= U1 

- Cl 

+ 

(G l 

SfMAH. 

1 

/ 4. 

* IUCS 

-Ul) 

VW2 

= U4 

+ C4 


(G i 

(MMA+1. 

) 

74. 

*(UCS 

-U4) 

IF 

(VW1 , 

.GE. 

U) 

GO 

TO 

5 





IF 

(UCS , 

• GEr 

UI 

GO 

TO 

10 





IF 

( VW2 . 

.GE. 

U) 

GO 

TO 

15 






P9 = P 4 
U8 = U4 
RHOB - PHO 4 
VB = V4 
RETURN 
PB = P1 
UB = U1 

P.H0B = RHD1 

VB = VI 
RETURN 
PB=PCS 
UB=UCS 

RHOB = RHOl* 1PCS/PH ** (1 ./GAMMA) 

VB = VI 

RETURN 

pb=pcs 

UR=UCS 

RHOB = RH04MPCS/P4 ) * * (1. /GAMMA) 

VB = V4 

RETURN 

END 
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APPENDIX B 


THE SHELL CODE 

As part of the evolution of GODUNOV, comparisons were made 
with another two-dimensional time-dependent hydrodynamic code called 
SHELL. SHELL was created several years ago for nuclear detonation 
calculations. It was intended originally in the present investigation 
to use SHELL in the focal region of the calculation. However, it is 
inappropriate for sonic boom studies because it tends to smear out 
weak shocks with relative overpressures less than 0.07. Nevertheless, 
SHELL was useful in checking out GODUNOV in certain test calculations 
(see Chapter 3) . 

A flow chart for SHELL is shown in Figure A-2. The input 
for SHELL is generated by an auxiliary program called CLAM. In the 
course of one time step in SHELL, the fluid properties at the center 
of each cell in the grid are updated in two phases. In the first 
phase, the conservation equations are solved with the convective terms 
neglected. In the second phase, material is allowed to flow across 
cell boundaries and transport mass, momentum, and energy. A detailed 
discussion of the calculation procedure can be found in General Atomic 
report number GAMD-5580, "OIL, A Continuous Two-Dimensional Eulerian 
Hydrodynamic Code," 1965, by W. E. Johnson. 
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FIGURE 1-2. GEOMETRIC ACOUSTICS DESCRIPTION OF FOCUSING, 
SHOWING RAYS AND CAUSTIC CUSP 


2 




7 ' 
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FIGURE 1-6. FOCUSING OF A SHOCK WAVE FROM A DIVE MANEUVER 




FIGURE 1-7. SONIC CUT-OFF FOR AN ACCELERATING AIRCRAFT 
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FIGURE 3-3. RIEMANN PROBLEM AT A CELL BOUNDARY 
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FIGURE 3-4. AXIAL SYMMETRY CONDITIONS IN GODUNOV 

TOP: VERTICAL CELL-BOUNDARIES AT AXIS 
BOTTOM: INCLINED CELL-BOUNDARIES AT AXIS 
(x DENOTES CENTER OF CELL) 
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40 80 120 160 200 240 280 

z (meters) 

FIGURE 3-5. FORE AND AFT SHOCK PROFILES OF AN N-WAVE PROPAGATING TO THE RIGHT 
(VERTICAL CELL-BOUNDARY SYMMETRY CONDITION) AT FIVE DIFFERENT TIMES 
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FIGURE 3-6. PRESSURE PROFILES IN AN N-WAVE (AP n /p n =0.1) (GODUNOV CALCULATION) 


1.2x10° t= .187 sec 

I (400 cycles 
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FIGURE 3-7. PRESSURE PROFILES IN AN N-WAVE (Ap_/p_ = 0.1) (SHELL CALCULATION) 






FIGURE 4-1. SHOCK SEGMENTS AND NODES IN WHITHAM 
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FIGURE 5-1. SHOCK PROFILES OF A POLYNOMIAL-FRONT N-WAVE (Ap/p = 10' J ) (ACTUAL FOCUS IS AT z=79 METERS) 





60 


70 


z (meters) 


80 


90 


FIGURE 5-3. 


PRESSURE PROFILES AT THREE DIFFERENT RADIAL POSITIONS 
JUST PRIOR TO FOCUSING FOR A POLYNOMIAL-FRONT N-WAVE 


(*p 0 /p 0 - i°' 3 ) 
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FIGURE 5-4, RELATIVE OVERPRESSURE VERSUS AXIAL POSITION OF FORE 

SHOCK FOR A POLYNOMIAL-FRONT N-WAVE (ad /p = 10' 3 ) 

o o 


90 







.013 
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FIGURE 5-7. PRESSURE PROFILES AT AXIS AT THREE DIFFERENT TIMES FOR A GAUSS IAN -FRONT N-WAVE. 
Up n /P n = 10' 3 , Ap/Ap n = 13.0) 



FIGURE 5-8. FORE SHOCK PROFILES OF A RELATIVELY STRONG GAUSSIAN-FRONT N-WAVE 
(Ap 0 /p 0 = 0.90; FOCUS OCCURS AT t = .037 SEC AT z = 38 METERS) 









(W) X 
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FIGURE 5-10. SINGLE RELATIVELY STRONG CONCAVE SHOCK (GODUNOV CALCULATION) 



fore shock 
aft shock 



97 . 


FIGURE 6-1. FORE AND AFT SHOCK PROFILES OF AN N-WAVE REFRACTED BY A COLD-SPOT AND 
SUBSEQUENTLY FOCUSED (r t = 50 METERS, TEMPERATURE RATIO = 0.5) 
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FIGURE 6-2. FORE AND AFT SHOCK PROFILES OF AN N-WAVE REFRACTED BY A COLD-SPOT AND 
SUBSEQUENTLY FOCUSED (r spot = 150 METERS, TEMPERATURE RATIO = 0.5) 
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FIGURE 6-3. RELATIVE OVERPRESSURE VERSUS AXIAL POSITION OF FORE SHOCK FOR N-WAVE 
REFRACTED BY A COLD-SPOT. ('Ap /p = .0167) 
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FIGURE A-2. SHELL FLOW CHART 















TABLE 5-1 


FOCUS FACTORS FOR POLYNOMIAL-FRONT 
N-WAVE WITH DIFFERENT GRIDS 


Number of Points 

Ap max 

radial x axial 

Ap o 

20 x 7 

6.4 

100 x 7 

12.4 

200 x 7 

13.1 

400 x 7 

12.9 

100 x 12 

16.2 

200 x 12 

18.8 

100 x 20 

18.7 

* 50 x 7 

13.1 




















TABLE 5-2 


FOCUS FACTORS FOR GAUSSIAN-FRONT 
N-WAVE WITH DIFFERENT GRIDS 



Number of points 
radial x axial 

A *“max 

A P 0 

20 x 7 

3.0 

* 50 x 7 

9.4 

* 50 x 20 

12.4 

* 50 x 50 

13.0 


1.05 (non-uniform spacing) 
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